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Chapter I. 


Project Summary 


1. Research Objectives 

The ultimate goal of the project is to theoretically investigate the effect of nonequihbnum. phase 
transitions on vapor and vapor-gas bubble dynamics in acoustic fields to determine the feasibility 
of the measurement of the accommodation coefficient using bubble dynamics in acoustic fields. The 
technical objectives include: 

1. Using a spherical model of a bubble in an isotropic acoustic field including thermal and 
diffusion effects in the liquid and in the vapor, surface tension, liquid viscosity and com- 
pressibility, evaluate the effect of nonequilibrium phase transitions on vapor and vapor-gas 
bubble dynamics . 

2. Evaluate the influence of the accommodation coefficient on the bubble translational motion 
in a standing acoustic waves and on rectified heat transfer. 

3. Conduct a parametric study to evaluate the mnge of parameters which can be used in design 
of an experimental setup for determination of the accommodation coefficient. 

2. Summary of Research 

The research effort has successfully achieved most the objectives set in the proposal: 

• We developed a comprehensive model of a spherical bubble in an acoustic field including the 
effects of nonequilibrium phase transitions, heat and mass transfer in the liquid and gaseous 
phases, surface tension, liquid viscosity, and compressibility. 

• We developed and implemented asymptotic and numerical methods for solution of the prob- 
lem and incorporated them into codes convenient for parametric studies. 

• We conducted parametric studies of vapor and vapor-gas bubble dynamics m acoustic fields 
and determined conditions in which the effect of nonequilibrium phase transitions on bubble 
dynamics is appreciable and can be used for measurements of the accommodation coefficient. 
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• We uncovered several physical effects which can be used for determination of the accommo- 
dation coefficient using bubble dynamics. These include, for example, low frequency bubble 
mean position/radius oscillations in standing acoustic waves, the possibility of stabilization 
of the mean bubble radius in acoustic fields, the existence of multiple threshold and stable 
equilibrium states of the mean vapor bubble radius in acoustic fields and their dependence 
on the accommodation coefficient, parameters of the acoustic field, and ambient conditions. 

• It is commonly thought that nonequilibrium phase transition affects bubble dynamics only 
in high-frequency fields (at least tens of kilohertz for water; some authors considered even 
the megahertz range). We, however, found that the range of bubble dynamics sensitivity 
depends on the bubble size, initial conditions, content of the inert component, and other 
parameters, which can be selected in experiments to provide measurements over a broad 
range of acoustic frequencies, and ambient conditions. We found that nonequilibrium phase 
transitions may affect the dynamics of vapor and vapor-gas bubbles in a broad range of 
acoustic frequencies and bubble sizes (for water at atmospheric pressure, 1-100 kHz and 
10-10000 fjb m, respectively). 

• We found a strong effect of initial conditions on rectified heat transfer for an initial stage 
that can span millions of cycles of bubble oscillations. 

• We found that gravity can substantially affect the bubble dynamics in acoustic fields. 
Regimes of bubble dynamics in standing waves substantially depend on the magnitude of 
gravity. Reduced gravity in general is beneficial for measurements of the accomodation 
coefficient. 

3. Problems and Future R&D Objectives 

The following research issues still need to be addressed prior to designing an instrument for 
measurement of the accommodation coefficient. 

1. Validation of the developed codes by comparison with available experimental data. 

2. Modification of the developed codes to include convective heat transfer due to bubble tians- 
lational motion in a standing wave. 

3. In the space of parameters, determination of the regions of spherical shape stability /applicability 
of the spherical bubble theory. 

4. Evaluation of the effect of bubble nonsphericity on the measurements of the accommodation 
coefficient. 

5. Determination of the optimal parameters/regimes for measurement of the accommodation 
coefficient. 
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4. Potential Applications 

The research seeks to develop a practical diagnostic technique for measurement of the accommo- 
dation coefficient of various substances in various conditions. The technique is based on measure- 
ments of bubble size and position in acoustic fields and subsequent processing of the results of 
measurements. Dependencies of the accommodation coefficient on temperature, surface contam- 
ination, and other parameters can be established using the proposed technique and utilized for 
diagnostic purposes and characterization of liquid/vapor interfaces. Measurement of the accom- 
modation coefficient is important for proper modeling of many natural and technological processes 
including boiling of liquid metals, film boiling, vacuum vaporization, explosions, aerosol mechanics, 
meteorology, and others. 
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Chapter II. 


Introduction 


Nonequilibrium liquid/vapor phase transformations occur in a wide variety of natural and tech- 
nological processes. In addition to the problem of interest here, these include evaporation and 
condensation of high velocity jets [8], film condensation [9], growth of small droplets in clouds 
[10, 11], sound propagation in vapor-droplet systems [12, 13], nonlinear dynamics of vapor bub- 
bles and condensed droplets in acoustic fields [14, 2], and laser vaporization [15]. 

Studies of processes with non-equilibrium phase transitions are of great practical importance, 
because in many advanced technologies it is necessary to predict and control material behavior 
under extreme conditions. However, the use of most theories and models is limited by the lack of 
reliable data on material properties, especially coefficients describing nonequilibrium vapor /liquid 
transformations. Such data can be obtained only from experiments. The experimental facilities 
used for these measurements are usually complex, expensive, and do not provide repeatable re- 
sults. For example, data reported on the accommodation coefficient for water obtained by various 
experimenters during this century vary from 6T0 -3 to 1, and experimental data for mercury aie 
in the range from 5T0” 4 to 1. 

The reason for such a wide range of results is that the accommodation coefficient is very 
sensitive to the conditions of the experiment. In order to address this shortcoming, the dependence 
of the accommodation coefficient on temperature, on the concentration of other species, and on 
other parameters near the interface should be determined. If a reliable and accurate accomodation 
coefficient measurement technique is available, such dependencies can be found experimentally and 
tabulated. This will provide a firm basis for modeling of nonequilibrium phase transitions. 

In recent microgravity experiments conducted by the European Space Agency using the Ger- 
man drop tower and the Bubble, Drop, and Particle Unit ESA multi-user facility for fluid physics 
experiments operating onboard IML2 (1994) and LMS (1996), an attempt to determine of the ac- 
commodation coefficient of refrigerants Rll and R113 was made by Picker and Straub [16] using 
observation of vapor bubble dynamics. Typical times of bubble radius variations in these exper- 
iments were of order 1 s. These times are several orders of magnitude larger than the interface 
temperature relaxation times and the evaporation/condensation in these experiments occurred in 
almost equilibrium conditions (no dependence on the accommodation coefficient). This explains 
the great dispersion of the values obtained for the accommodation coefficient (from 0.9*10 to 0.7 
for Rll and from 8*10“ 3 to 1 for R113 [16]). 
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It is known that during each cycle of bubble oscillation there are two phase transition stages, 
evaporation and condensation. For very small amplitude oscillations the amounts of evaporated 
and condensed liquid are approximately the same, and the vapor bubble oscillates about an equi- 
librium value. However, at larger amplitudes the nonlinearity of the bubble dynamics causes a 
difference between the amounts of evaporated and condensed liquid. This difference is small duiing 
one period of oscillation, but leads to slow-timescale dynamics of the average bubble size. This 
is called “rectified heat transfer”. The same effect for gas bubbles growing due to mass diffusion 
is known in literature as “rectified diffusion”. For vapor/gas bubbles it can be called “rectified 
heat and mass transfer”. In some regimes the average bubble size can reach an equilibrium value, 
and the bubble can experience stable oscillations. In this state the amounts of evaporated and 
condensed liquid over the period are equal even when taking into account all nonlinear effects. The 
equilibrium mean radius can be defined as the mean radius separating ranges of bubble growth 
and shrinkage. The mean equilibrium radius can be stable (in this case bubbles of sizes slightly 
larger than the equilibrium radius shrink and bubbles slightly smaller than the equilibrium radius 
grow) or unstable (opposite situation). In a standing acoustic wave the bubble also experience a 
slow drift under the action of the primary Bjerknes force, which is also a nonlinear effect. The 
present study shows that the effects of bubble drift and rectified heat transfer are strongly coupled. 

Accounting for these nonlinear effects in a measurement technique has additional advantages. 
First, these effects depend strongly on the value of the accommodation coefficient. Second, the 
characteristic times of the growth of the bubble average radius and of the bubble drift m weak 
acoustic fields are much larger than the bubble period of oscillation. Also the corresponding 
spatial scales (of the order of the bubble radius and the acoustic wavelength) are much larger than 
the amplitude of the bubble radius oscillations. This makes these measurement much easier to 

accomplish. . 

Successful completion of this effort will enable the refined theory and codes to be used tor 
development of a simple system for measuring the accommodation coefficient. Such a system 
could be used not only in fundamental studies of kinetics of phase transitions, but also could have 
good practical applications, for example for detection of extremely small amounts of contaminants 
present in pure substances, because of the high sensitivity of the accommodation coefficient to 
the contaminants. To establish such procedures the dependence of the accommodation coefficient 
on various contaminants and ambient conditions should be determined first, which also requires a 
good measurement technique for this coefficient. There could be substantial benefits from using the 
accommodation coefficient measurement technique for controlling production of pure substances 
such as semi-conductors, high purity chemicals and isotopes. This technique can be modified or 
the case of arbitrary aggregate states of matter, and allows consideration of gas-solid, liquid-liquid 
and liquid-solid non-equilibrium phase transitions (e.g. the crystal growth from melts in acoustic 
fields under the effect of rectified heat transfer as was considered in [17]). 
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1. Historical review of measurements of the accommoda- 
tion coefficient 

The kinetic model of evaporation and condensation was first proposed by Hertz [18] and Knudsen 

[19] . If the difference between the saturation temperature T s and the temperature of the interface 
T a is not too large, the rate of evaporation £ can be written in the form of the Hertz-Knudsen- 
Langmuir equation. The accommodation coefficient (3 is also known as the evaporation or the 
condensation coefficient. Some authors discriminate between the condensation and accommoda- 
tion coefficient, while other authors use the same value for both terms. /? is a dimensionless 
thermophysical parameter ranging from 0 to 1 which is a property of the liquid-gas interface (such 
as the surface tension coefficient). 

The history of accommodation coefficient measurements shows many problems, as can be 
seen in the publications cited in this section and the fundamental monograph of Volmer [1]. All 
these works and the experimental data show a wide range of possible values of (3 and show the 
importance of considering the physicochemical hydrodynamics near surfaces where nonequilibrium 
condensation occurs. The following examples illustrate this point. 

The first experiments of Hertz with mercury [18] gave the value (3 = 0.11. Knudsen [19] also 
experimented with mercury, and his first experiments gave the value (3 = 0.0005. He explained 
this low value of the accommodation coefficient by surface contamination. After improving the 
conditions of the same experiment he obtained (3 = 0.11. Finally, experimenting with falling 
mercury drops with continuously renewing surfaces he found (3 = 1. 

Numerous experiments of Langmuir and coworkers performed between 1913 and 1934, showed 
that the accommodation coefficient for high temperature boiling metals such as wolfram, molyb- 
denum, platinum, nickel, iron, copper, and silver is close to 1. Further analysis of these and other 
data showed that the value of (3 for all of the mentioned metals is between 0.25 and 0.33, except 
for platinum, for which (3 = 1 could be true. 

The accommodation coefficient for water was measured in the studies of Alty and coworkers 

[20] . To determine the interface temperature they used a technique based on the measurement 
of the surface tension and found (3 = 0.036. In [9] the water accommodation coefficient was 
estimated from measurements of filmwise condensation of steam at low pressure on a vertical flat 
plate and (3 was found between 0.45 and 1. On the other hand, recent measurements of the water 
accommodation coefficient using a cooled wall expansion chamber [11] showed values of (3 that 
decreased from 1 for 1 /mi droplets to 0.006 for 15 fim droplets. This, perhaps, could be due 
to some diffusion of trace contaminants or inert components to the surface. Nevertheless, the 
experimenters that took precautions to ensure system and sample purity were not able to identify 
trace contaminants. The authors of [21] also found a difference in the water accommodation 
coefficient between a fresh surface (measured (3 is 0.2) and a “stagnant” surface (measured [3 is 
O.O 38 ). 

In Fig. II-l we present some reported data on accommodation coefficient measurements. We 
plotted the ratio of the maximum reported to the minimum reported value of the accommodation 
coefficient for five different substances. It is seen that these ratios can be of order 10 2 or even 10 . 
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Figure II-l: The ratio of the maximum to the minimum measured values of the accommodation 
coefficient for five different substances. 

2. Review of publications related to vapor bubble dynam- 
ics 

The acoustics of vapor bubbles and acoustic vapor cavitation were intensively studied theoretically 
and experimentally in the 1960’s and 1970’s. Many of these studies were related to the design of 
cryogenic bubble chambers for registration of the tracks of charged particles [22] and measurements 
of the tensile strength of liquids [26]. More recent applications include acoustic enhancement of 
boiling in microgravity [27, 28] and the use of bubble dynamics for determination of liquid- vapor 
interface properties [4]. 

Wang [29], Khabeev [30], and Fanelli et al [31] performed linear analyses of forced vapor and 
vapor-gas bubble oscillations and showed a strong difference in acoustic properties of vapor bubbles 
and bubbles of non-condensable gas. In addition to the primary resonance, vapor bubbles exhibit 
a second resonance in acoustic fields corresponding to smaller sizes for a given frequency, which is 
known as the condensation-evaporation resonance. This resonance was first reported by Finch and 
Neppiras [32]. Hsieh [33], Marston [34], and Hao and Prosperetti [35] provided physical insight 
into the second resonance. 

The effect of rectified heat transfer on vapor bubbles was investigated theoretically by several 
researchers [3, 5, 35, 36, 37, 39]. Marston and Greene [40] observed stable oscillations of bubbles 
for several seconds in liquid helium-I. Our recent study sponsored by NASA shows that there can 
exist multiple threshold and stable equilibrium mean radii of vapor bubbles m isotropic acoustic 
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waves [5]. Hao and Prosperetti [35] and our study [5] have determined that the development of a 
“slow” thermal boundary layer in the liquid at large time is a significant mechanism of rectified 
heat transfer. This means that convective heat transfer can be important for bubbles moving in the 
liquid. Theories and computational schemes for modeling such processes are available [41, 42, 43]. 

There are very few publications available on experimental studies of rectified heat transfer. For 
example, Akulichev et al [17] mention some experiments in cryogenic liquids without providing 
quantitative information. A very recent publication of Ohsaka and Trinh [44] reports results on 
water vapor bubble growth rate measurements, which perhaps is the first publication of such 
data. This paper shows that experiments with vapor bubbles in acoustic fields are challenging. 
Particularly, experiments for “vapor” bubbles were performed at atmospheric pressure and at a 
temperature of 80°C. Such bubbles consisted approximately of 40% water and 60% air (by weight). 
Our computations show that such bubbles differ from both pure air and pure vapor bubbles, and 
a theory of rectified heat and mass transfer to vapor-gas bubbles should be applied in this case. 

The above mentioned theoretical studies deal with spherical bubbles. However, due to the 
parametric resonances between volume and shape modes, a stable spherical bubble shape can be 
realized only at small amplitudes [45, 46, 47]. Other effects, such as acoustic streaming, translatory 
bubble motion, and gravity, can also influence the results [48]. To describe shapes of levitating 
bubbles and drops in standing acoustic waves variational and other methods were applied (e.g. 

[49, 50]). 

Forces acting on the bubble in acoustic fields were studied by several researchers. The major 
force acting on a bubble in an oscillating pressure field is the primary Bjerknes force [51], which 
is the time average over a period of the product of the bubble volume and the acoustic pressure 
gradient. This force is proportional to the energy of the acoustic field. In normal gravity conditions, 
depending on bubble size, frequency, and amplitude of sound, it can exceed the gravity force [52] 
and the bubble can levitate [53] or be positioned in the center of a spherical flask such as in 
sonoluminescence experiments [54]. Other important forces include the added mass force and 
viscous drag force. 

Note that the value and sign of the primary Bjerknes force depend on the response of the 
bubble volume to the acoustic excitation. Since there exist a substantial difference in resonance 
properties of gas, vapor, and vapor/gas bubbles, the primary Bjerknes force is different for these 
three cases. At higher frequencies it depends on the kinetics of phase transition due to its influence 
on the bubble resonance [6]. Thus the model predicting acoustic forces on the bubble and the 
resulting bubble motion and shape deformations should include an accurate consideration of bubble 
forced oscillation including heat and mass transfer inside and outside the bubble, kinetics of phase 
transitions, liquid inertia, compressibility, viscosity, and surface tension. 


9 




Chapter III. 


Statement of the Problem 


We will consider the dynamics in acoustic fields of a bubble filled with the vapor of the host 
liquid or a mixture of the vapor and an inert gas. We consider that the boiling point of the 
inert gas corresponds to substantially lower temperature than the liquid temperature at the same 
pressure. A mixture of water vapor and air is an example of a such system. We also consider two 
configurations of the acoustic field: a) an isotropic field; b) a standing wave. 

1. Model of Vapor-Gas Bubble in Isotropic Acoustic Field 

Consider a spherically-symmetric model of a vapor-gas bubble in an isotropic pressure field, with 
the wavelength much larger than the bubble size, uja <sjC C, where u> is the circular frequency, a is 
the bubble radius, and C is the speed of sound in the liquid. For a viscous liquid and inviscid gas 
the mass, momentum, and energy conservation equations at the interface can be written in the 
form [38, 58, 59]: 

Pl ( a — U>ia) = pga (d — W ga ) = f , f = ft, + ft, (IH'l) 

n£ = -Pg + f faga ~ Wla) + 2~, (HI-2) 

1 O • 2 O' CL /TTT 0\ 

n^Wla — qia + —£w* a = — + fit + d" H ( ) 

Here p,w, and q are the density, radial velocity, and heat flux, p and n rr are the pressure and 
radial component of the stress tensor, and f,cr, and l are the rate of phase transition from liquid 
to gas, surface tension, and heat of phase transformation. Subscripts l and g refer to liquid and 
gas, respectively, and subscript a denotes parameters on the interface. We assume that the gas 
consists of two components, vapor and inert gas, which parameters are marked with subscripts v 
and i. 

The dynamic equation describing forced radial oscillation of a bubble of variable mass bubble 
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in slightly compressible liquid can be found in [60, 5]: 



aw la + 2 ( 

a a d\ 
1 + C + Cdt) 


1 Wla\ - 1 2 

1 - aw«, - -w la 

) II£ + Poo (t) + ~IMWia . 
i a 


(III-4) 


Here pi is the liquid viscosity and p^ (i) is the forcing pressure. 

In the general model we accept that the gaseous phase is a perfect gas mixture described by 
the following thermodynamic relations: 


Pi — RTg Inpi , p v — RTg]np V: fig — pi p v . (III-5) 

Pi == PiRiRg , Pv = PvRyRg > Pg Pv 4" Pi i (III-6) 

Ci = — , C v = —, pg=pv+pi, Ci + Cv = 1, (HI-7) 

iii = H/M, R V = R/M V , R g = c i R i + c v R v . (HI-8) 


Here p is the chemical potential, T the temperature, c the mass concentration, R the gas constant, 
and M the molecular weight. 

The resulting mass fluxes of the components from the liquid to gaseous phase, £, a.t the constant 
interface temperature T a can be described by [1]: 


6 = ^ 9) ~^ l) = 




A 


{l ^ 9) l4 (Ta) 

exp 


V^wRjTa 

$-* 9) 


y/2 ^Rjr a exp RT a V27T 7WT C 


RT a y/2n RjTa 

Pi (To) 


/* (T a ) 

“'P-RTr’ 

pu (To) 


exp ■ 


RTa 


where P are the proportion coefficients and the superscripts near j and 0 denote the direction of 
the mass flux. At thermodynamic equilibrium, j g = j v = 0, p v = p v3 ( T a ) , and p, = C[ a H. Thus 
= pi and Pv~* 9 ^ = Pv 9 ~* 1 ^ = /?„. Using the above expressions for the fluxes and 
chemical potentials we have 


(IIM) 

The last relation is the well-known Hertz-Knudsen-Langmuir equation describing non-equilibrium 
evaporation of a one-component liquid [1]. The coefficient P v is the vapor accommodation (con- 
densation) coefficient and can be treated as the fraction of vapor molecules hitting the interface 
which condense. By analogy Pi can be called ‘inert gas accommodation coefficient’ and can be 
measured from experiments on nonequilibrium dissolving. 

From a linear analysis described in Chapter 4 and in Chapter 7 we found that Henry’s law, 

Pia = c la H, (HI-10) 

is applicable for description of bubble dynamics over a broad range of frequencies, while for the 
vapor component the nonequilibrium evaporation/condensation is important. 
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For spatially uniform pressure the diffusion flux of the inert component, j g , and the heat. flux, 
q g , can be expressed according to [56]: 



k r dT g \ 
T, dr)' 


dT 

Qg = ^ ^cCpgTgjg, 


(III-ll) 


k C = 


1 

CpgTg 





k T R 

CiC v Cj 


P9 


(III-12) 


where D and A are the mass diffusivity and thermal conductivity, kx is the thermal-diffusion 
ratio, and C pg is the gas specific heat at constant pressure and concentration. Expression (III-12) 
is obtained using the above model of the gaseous phase (III-5)- (III-8) , where fi g is considered as 
a function of p g ,T g , and c*. Note that in limiting case of one-component gas, QCy —* 0, we have, 
kx —* 0, while k c remains to be a finite quantity. The quantities k c or kx can be found from 
corresponding tables or can be evaluated using the formulas following from the above definitions 
and found in Ref. [61]: 


a Mj R _ 105 Mj - M, . 

MiCv 4- M v Ci Cpg ’ a 118 Mi + M v T 


(III-13) 


where a is the thermal-diffusion constant, and k>p has been evaluated for several intermolecular 
force-models. For rigid elastic spheres k^. = 1 and we took this value for our computations. 
Details for computation of the thermal-diffusion ratio based on the first, second, and third order 
gas-kinetic theories can be found in Ref. [62]. 


Now we can represent the mass and energy conservation equations in the form: 

dp g 13m \ _ / dci dci\ Id , 2 >. 

„ (dTg 8T g \ ... 1 d , 2 . 


(III-14) 


In these equations we assume that the total gas pressure is spatially uniform, which is justifiable 
when the velocity of the bubble wall is much smaller than the speed of sound in the gas. Note 
that at the same time the partial pressures of the components depend on the radial coordinate 
due to dependence of the inert gas concentration on this coordinate. 

Assuming that the mass concentration of the dissolved inert gas in the liquid is small, Q <?C 1, 
we can neglect the effect of thermal diffusion in the liquid. The effect of barodiffusion is also 
negligible, since the liquid is almost incompressible. Therefore, we can represent the mass and 
energy conservation equations in the liquid in the form: 


f dci a 2 Wi g dci \ 

P\dt+ r 2 dr ) 

dTj a?wia dTj \ 

dt r 2 dr ) 



1 d , 2 ,\ 

r 2 Qr ( r 3 l ) ’ 

1 9 ,2 v 12 fnwf a a 4 

( r 9/) + 3 


n d ° l 


Qi = 


. m 

~ Xl ~RZ' 

or 
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(III-15) 

(III-16) 


where Ci is the liquid specific heat. 

Assuming that the temperature jump at the interface is negligibly small (which is true for not 
too high rates of evaporation or condensation), the following boundary conditions can be imposed: 

W g\r=a = ^pl r =a = ^lr=a = T[\ r=00 = Tqo, (III-17) 

C ilr=a ~ C /| r=a = Cla > C l r =oo = C °°* (III-18) 

The diffusion mass fluxes at the interface can be represented as 

jga = jla ~ £t- 


To specify p V3 {T) the Clausius-Clapeyron equation can be used: 

dpvs _ iv f i _ IV 1 
dT ~ T \p vs (T) pj • 


(III-19) 


where p vs is the vapor density on the saturation line. In the present study for simplicity we neglect 
the dependence of a,H,Xg,Xi,D g , and Di on the temperature. 

Note that if the heat and mass fluxes are known, then these equations together with the 
Rayleigh-Plesset equation (HI-4) form a closed system. These can be found by solving corre- 
sponding problems of heat and mass convective diffusion with boundary conditions Ci a and C; a for 
concentrations and T a for temperatures inside and outside the bubble. 

We also need to specify the function ( t ) , which for an acoustic field of amplitude Pa and 

circular frequency u can be written in the form 


Poo ( t ) = Pooo 4- Ke{P A e wt }. 


(III-20) 


2. Model of Vapor Bubble in Isotropic Acoustic Field 

The case of pure vapor bubbles is a limiting case of the vapor-gas bubbles at small concentrations 
of the inert gas. However, it is important to consider this case in parallel with the case of vapor-gas 
bubbles, due to substantial model simplifications that can be obtained for one-component bubbles. 
This limiting case can be used for verification of the general results for two-component systems. 
Simplifications for pure vapor bubbles compared to vapor-gas bubbles include: 

• Simplification of kinetics of phase transitions; 

• Simplification of bubble thermodynamics. Availability of a simplified energy integral for 
vapor; 

• Absence of mass diffusion; 

• Simplification of boundary conditions. 
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Equations describing the vapor bubble dynamics can be represented in the form: 


pi(a-wi a ) = p va (a- w va ) = £, 

n" = -p v + ( {u>va - ™la) + 2- , 

1 1 2 2crd 

IIZw la - qia + ^wf a = -p v w va - q va + ^ w va + ^ + O' + — 

(i - aw Ja + 2 (l - aw ia - 2 W ia 


1 / a a d \ 
= ~J l { 1 + C + Cdt) 


+Poo (0 + -/^«7u. 


a 


Pd — Pv-^v^gi 


(III-21) 

(III-22) 

(III-23) 


(III-24) 


^ (Ta) P ” J ' 


(III-25) 




= 0, 


1 ^ , 2 x 

~7d? (r ?v) ’ ? ’ = 


(III-26) 


^ ;m , a 2 w/ a <97} 


) = -^1 (’■ 2 *) + 


2 „4 


12 pi wf a a‘ 

r 6 


. m 

q i- 


(III-27) 


^U = «va, T,u = T,u* r .. r^-Too. (HI-28) 

Since for one-component bubble content the pressure of the vapor is spatially uniform (for a 
two component system the inert gas and the vapor pressure are not uniform due to gradients of 
concentration) the equations for vapor inside the bubble have the following integrals [59, 57, 63]. 


_ (7« ~ l)*v&Tv _ r Pv (HI-29) 

Wv 7 v pv dr 3 t v pv ' 

ap v + 3 'YvPvWva + 3(7,, - 1)9™ = 0, (III-30) 

where 7„ is the ratio of the vapor specific heats. These integrals simplify solution of the problem 
since they explicitly express the vapor velocity through the temperature gradient and vapor pres- 
sure, and connect variation of the pressure in the vapor with integral parameters characterizing 
the bubble (bubble radius, radial bubble wall velocity, and the heat flux through the surface). 
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3. Model of Small Bubble in Standing Acoustic Wave 


Consider the pressure and velocity fields, p (x, t) and Uoo (x, t), of a plane standing acoustic wave 
in a liquid: 


Poo («, t) 
Uoo (x, t) 


Pooo + Pa cos u)t sin kx, 


Pa 

PiC 


smut cos kx, 



where k is the wavenumber. We limit our analysis to bubble radii small compared to the wave- 
length, ka <C 1. The acoustic field in the liquid described by Eq.(III-31) is obtained using linearized 
equations for the liquid motion. This approximation assumes that 


Uo 


dU 0 


dx 


< 


dU a 


dt 


(III-31) 


or Pa PiC 2 , which holds for a wide variety of situations. 

The bubble position, Xb(t), can be found by solving the following equations for bubble motion: 



dxb 

dt 


--ira 3 (VpJ,„, - p,g) - m b g + 

d(a?U x U,) i(a>U b ) 
dt dt 


r pi 


+47rhT M pia (Uoo ~ Ub ) , 


+ 


U b , 


drrib 

dt 


= 47ra 2 £. 


Here Ub and mb are the bubble velocity and mass, £ is the rate of evaporation, g is the gravity 
acceleration, and K ^ is the viscous drag coefficient, which depends on the Reynolds number of 
the relative bubble motion. According to the Levich formula at high Reynolds numbers = 3. 
For small Reynolds numbers it can be set to = 1 or = 1.5 for a liquid with and without 
surfactants. Equations (III-32) are projections of three dimensional equations of motion on the 
direction of bubble motion. To apply one-dimensional equations we assume that the vector of 
gravity acceleration is co-linear with the direction of the acoustic wave vector and the direction of 
the bubble motion. 

We also modify the dynamic equation of bubble motion as 


(l ~ o.w la + 2 awi a - w 2 la = 


(III-32) 


a a d\ 

\ 1 + c + c7t) 


HZ + Poo + AfliWia/ a - Jpi (U b - UcoY 


to include variation of the liquid pressure at the bubble location and the bubble relative motion. 
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4. Vapor-Gas Bubble Equilibrium in the Absence of an 
Acoustic Field 

Vapor, gas, and vapor-gas bubbles can grow or shrink in liquids without any acoustic field. This 
happens due to phase transformations such as condensation or dissolution of a thermodynamically 
nonequilibrium bubble. A bubble in the bulk of the host liquid also cannot be in equilibrium under 
action of the gravity force. The buoyancy force will cause bubble motion in the absence of an 
acoustic field. If the thermodynamical and gravity driving forces are sufficiently high then the 
external acoustic field may not affect the bubble dynamics and this case can be studied separately 
without complication of the problem by acoustic action. In the present study we consider the case 
when the effects of an acoustic field are comparable to or larger than the effects driving the bub e 
dynamics in the absence of the acoustic field. We limit our study to small amplitude acoustic 

fields, 


£*. = £ « i. <m- 33 ) 

V ooO 

Therefore in this case the effects of bubble thermodynamic and gravitational instability m the 
liquid are small, and we assume first that the bubble is in a state close to thermodynamic equi- 

For vapor bubbles this means that the mean liquid pressure is close to the saturation pressure 
at the liquid temperature far from the bubble: 


P «o= P ,ey + Ap, A = ^«l. (in-34) 

where A p is the liquid supercompression. For superheated liquids we have Ap < 0 1 and for 
subcooled liquids we have A p > 0. We allow Ap to be non-zero, but it should be reasonably small 
(we provide mathematical definition in the sections dedicated to the method of solution). Even if 
Ap = 0, the vapor bubble is unstable due to capillary effects. Indeed the pressure inside a bub e 
of radius a e in thermodynamic equilibrium is 


+ Ap+f a“- 35 > 

LLe 

The temperature of an equilbrium bubble should be the same as the temperature of the liquid, 
T x , and at the same time it should be the saturation temperature at the given vapor pressure 
(III-35). This means that the equilibrium bubble radius is 




2a 

Ap’ 


(III-36) 


and equilibrium bubbles can exist only in superheated liquids. Note that this is an unstable 
equilibrium. In this and other cases a vapor bubble will grow or shrink in the absence oi the 
acoustic field. We require that 
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< 1, 


(III-37) 



A 

Ap-I 

a 


to have comparable effects of the acoustic pressure variations. 

The pressure inside an equilibrium vapor-gas bubble is 

Pg = Pv + Pi = PooO + — > PooO = PsiTeo) + A p, (III-38) 

CL e 

where the difference between the mean liquid pressure and the saturation pressure may be not 
small. In an equilibrium state the vapor-gas bubble has the same temperature as the liquid, 
T g = Too- If effects of inert gas dissolution are neglected, then the inert gas concentration in the 
bubble can be arbitrary, d 0 . In the equilibrium state we have p v = p^Tx,). Using equations of state 
(III-5)-(III-8) we can determine the equilibrium pressures of the inert gas and the gas mixture: 


Cio-R* irn \ _ CyoRy (q-< \ 

Pie ~ c v oR v Pa ^ r °°^' Pae ~ dioRv Pa{ 00)1 
Then we can determine from (III-38) the equilibrium radius: 


d>0 — I do* 


(III-39) 


2 ex 2u CyQ 

Pie A p doRiPs (Too ) f^yAp 


(III-40) 


If we take into account solubihty of the inert gas then its concentration inside an equilibrium 
bubble is not arbitrary, but is determined by the equilibrium relation between the inert gas pressure 
and the dissolved gas concentration in the liquid, c 0 0 (in our model, by Henry’s law): 


Pie — Cqo-H, de 


RvPi, 


RrPs {Too) "b R-vPie 


Pge — Ps (Too) *b Pit 


dp 


2(7 

Pie ~ Ap' 


(III-41) 


In any case a bubble of radius a ^ a e will grow or collapse in the absence of an acoustic field. The 
driving thermodynamic force is proportional to the difference between the actual and equilibrium 
gas pressures, and conditions (III-34) and (III-37) for vapor-gas bubbles become. 


A = l Pge Pooo| <l, <5 = — 

PooO PooO 

where relations (III-39) or (III-41) can be used to determine p ge . 

A parameter responsible for bubble instability due to the buoyancy force can be found by the 
following reasoning. In a standing wave the hydrostatic pressure drop, A p h , along the wavelength 
should be much smaller than the ambient pressure. Otherwise a small amplitude acoustic pressure 
cannot balance the hydrostatic pressure drop, and bubble motion will be not influenced by the 
acoustic field. This condition can be written in the form: 


2(7 

Pge PooO "h 


< 1, 


(III-42) 


A B= ftCs« li ( Aph „ ML „e!£l). (III-43) 

PooO OJP00O V U / 
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5. Length Scales and Dimensionless Equations 

The following characteristic length scales can be introduced for the problem considered. 


Lv= l -{- 

P W V Pi 


= ( — 

\P&> 


L a = 


Pg * 
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u 


(III-44) 




U = 


PI 



(1 - p)u>-J2'kR v T 00 ' “ u \ Hi 

where k, = A,/(p,a) and K ,. = are thermal diffusivities of the liquid and the 

gas, and stars denote characteristic values of parameters used for scaling. Lr g and -L X / a e 
characteristic lengths of temperature penetration into the gas and liquid, Lp is the characteristic 
primary resonance length, is the characteristic thickness of the viscous boundary layer m the 
liquid, L a is the characteristic capillary length, L c is the inverse wavenumber, and L 0 and L d 
are the characteristic lengths connected with the non-equilibrium phase transitions and viscous 

dissipation in the liquid. 

The following dimensionless parameters also can be introduced: 
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Note that the length scales written above in combination with the dimensionless parameters pro- 
duce new length scales having physical meaning. For example, characteristic lengths 



Lrg , = ( £i\ 1/2 _ J^ti_ 

JZT g ' cl W 


(III-45) 


represent the characteristic thicknesses of the mass diffusion boundary layers m g ^ 

liquid. The ratio of the bubble radius to any characteristic length scale listed in (I ) P 
8 independent dimensionless parameters. Taking into account that all but two of a {llbl3),e 
(HI-33), 6, (III-42), and *,(111-43) and parameters (IH-45) are also independent we can see that 
the dynamics of a vapor-gas bubble is controlled by 26 basic dimensionless parameters (!). For 
pure vapor bubble this number reduces to 15. 
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Chapter IV. 

Multiscale Technique For Bubbles in 
Isotropic Fields 


To obtain the equation for rectified heat transfer we use a multiscale asymptotic technique foi 
weakly non-linear oscillations of drops and bubbles [14, 2]. However, the method must be modified, 
since the effect of surface tension neglected in [14, 2] changes the rank of the system matrix for 
zero-mode of oscillation, and a straight-forward application of the technique is impossible. To 
demonstrate the method of solution we will start with a simplified system for pure vapor bubbles. 
Then we show how to include effects of bubble drift in standing waves and two-component effects. 


1. Transformation of Variables 


First, we reduce the number of variables by eliminating p v , w v, Pvai w vai an< ^ from the system. 
The governing equations can be represented in the following form: 

Pid - piWia - f = o, (rv-l) 
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We then transform — » (77, t) with p = r/a(t) to fix the moving boundary in the heat 

transfer problems: 


PvTooa 2 l dT v (7 - 1 )pyTy fa p v A dT v 

Ps {T^) Kv\_dt 7 Pv \ a 3 t p v ) 77 dp 


T v d / 2 dTv 
1 f dp dp 
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(IV-5) 
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a 977 


> ?=i 


(IV- 6 ) 


where /c* and /c v are the thermal diffusivities of the liquid and vapor (at ambient conditions). 

Next, we introduce dimensionless fast and slow time scales, to = an d £n = € n cjt, n = 1 , 2 ,... 
All unknowns are considered now as functions of this set of times and the temporal derivatives 
are represented as series: 


d ( d d 2 d \ 
dt=“Wo +€ m, +( dT, + --j- 

Finally, we expand the unknowns in the following asymptotic series: 


a(t) = a,o(ti,t2 , -.) [1 + ca,i {to,ti, ... ) 4 - ...] , 

w la(t) = WOo(ti,t 2 , ... ) [€U7i , 

Pv{t ) = Ps{Too) [1 +CPl (to, ill— ) + — ] 1 

T a (t) = Too [1 + zT\ (to, h , , 

£(t) = a;p t;oo ao(ti,t2, —) [ € Ci (^0,^1, •••) + —] > 

Qla{t) = a, roo [a (to, <!,.••)] 1 [ c 9i (to,tl,— ) + — ] > 

9 ua(t) = A^Toola^ti,...)] 1 [eri (to, ti,... ) + •••] i 

Ti(p,t ) = Too [l+€t*i (77, to, ti, , 

T v (p,t) = Tooil+tV! (77, t 0 , ti,.. .) + •..], 

where p voo is the vapor density at T v = Too and p v = p s {Too). 


(IV-7) 


(IV- 8 ) 

(IV-9) 

(IV-10) 

(IV-11) 

(IV-12) 

(IV-13) 

(IV-14) 

(IV-15) 
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2. Solution of Thermal Problems 


Formally, as e — > 0, the convective and source terms in the energy equations are small and, in the 
777 -th order approximation, we have the following inhomogeneous linear equations to determine 
the temperatures outside and inside the bubble, respectively: 


Qp dUm _ 1 d ( 2 dUm \ 

dto if dr] V dr] J 
_^_dvm _ Qp dp m _ l_d_ f i dv m \ 
L 2 Tv dt 0 L 2 Tv dt 0 V 2 dr] V dr l J 


= frr 


— Sm> 


(IV-16) 

(IV-17) 


Vml^o < 00 > (IV-18) 

where and g m (ri,to,t u -) are functions that depend on approximations of order 

less than m. 

Let us evaluate the range of frequencies, where this scheme is valid. Since for vapor bubbles 
the heat transfer in the liquid plays the major role, consider as an example the energy equation 
for the liquid (IV-6), for which the terms can be evaluated as follows: 


These are subject to the boundary conditions: 

= T m , Um\ r j =00 — 0) u mlj}=l T m , 
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piKiT] 2 dr] KiUdi 


(IV-19), 


Here Aa and AT are the characteristic variations of the bubble size and liquid temperature, di 
is the thermal boundary layer thickness, and U is the characteristic time. Consider fast bubble 
oscillations. In this case we have Aa ~ eooi AT ~ e ^ooj di ~ An A ~ u> 1 . At e ~ 0.1 the 
latter term in (IV-19) is small compared with the first term for w « = P/C/T*,/ Pi- Normally this 

limitation is not restrictive (e.g. for water and helium at atmospheric pressures we have u d /2it ~ 
10 u Hz). The convective term is related to the liquid motion generated by the moving bubble 
surface (the third term) and condenstation or evaporation (the fourth term). The third term (IV- 
19) is small at small e even at high frequencies, since in this case high temperature gradients of 
order eT^ j L ti are realized in a thin boundary layer r] — 1 ~ An/oo and (77 — 77 2 ) dTi/dr] ~ eT*,. 
The fourth term (IV-19) is small compared to the first term (IV-19) if epao < An- For e ~ 0.1,- 
p ~ 10~ 3 , and bubbles of typical radius 00 < 1 mm, the thickness of the thermal boundary layer 
in the liquid should be much larger than 0.1 pra. For water at atmospheric pressures this limits 
the theory to frequencies oj/ 27 t < 1 MHz. Note that the acoustic wavelength is of order 1 mm at 
frequencies o;/27r ~ 1 MHz, and the theory is limited by such frequencies for 1 mm bubbles. 
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We assume that all functions in (IV-16) 

^771 0 ?, ^ 0 , ^1 , •**) 

Pm (^0, tj , •••) 

where u mn and p mn are the complex amplitudes (the number of modes increases with the number 
of approximations because of non-linear generation of subharmonics). 

We should also notice that the first two orders of approximation of the present theory do not 
depend on initial conditions. This assumes that the bubble resides in the liquid for a long time 
and a quasi-steady slowly changing average temperature profile is developed near the bubble. In 
the m-th order approximation equation (IV-16) written for zero-mode 


- (IV-18) are periodic with respect to t 0 . For example 


Re \ ^ ] ^mn 0?> ^1 j ^2> •••) ® 


into 




n — 0 
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Re \ ^ ^ Pmn (^1,^2, •••) ^ 
v 71=0 


into 


(IV-20) 

(IV-21) 


1 d f 2 ^^77l0 'N 

rf dr] \ dr] J 
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(IV-22) 


shows that the non-stationary, x%f , and convective, Pe Q (tj - r\~ 2 ) terms are considered to 
be small compared to the conductive term ^ (V^f) and are passed to the approximations 
of order higher than m. Here Pe 0 = is the Peclet number of the slow bubble motion, and 

^ (we show later that nothing in the present theory depends on t\, so k ^ 2). This is correct 

from the point of view of the current formal asymptotic procedure, however, is questionable from 
the physical view point, because actual values of x an< ^ -^ e o may be not small. Bringing these 
terms into the m-th order approximation creates a mathematical problem, which solution currently 
is not available. Note that Fyrillas & Szeri [64] obtained a solution for transient slowly evolving 
fields considering the non-stationary term to be of the same order as the conductive term, but for 
constant a 0 (Pe 0 = 0). If significant growth or shrinkage of the bubble occur m some time scale 
t k (so dao/dt k ~ a 0 ), then in that scale Pe 0 ~ X, and the convective term is of the same order of 
magnitude as the non-stationary term. 

The assumption on periodicity of functions in the fast time scale is not valid for initial value 
problems at times of order of one period of oscillation. Fyrillas & Szeri [64] evaluated that time as 
a few periods of oscillation. The depth of penetration of temperature perturbations into the vapor 
and liquid due to oscillations of the bubble surface temperature is much smaller than the thickness 
of the thermal boundary layer evolving in slow time scales. That is why the stage of establishment 
of a periodical solution is not important when we consider slowly changing temperature fields. 
Since slowly changing fields depend only on slow variables the initial conditions for their evolution 
should also be formulated in slow time scales. In the present theory the initial conditions appear m 
the third-order approximation, which also requires a spatial matching procedure and is described 

in a separate section. _ , 

For each mode, Equations (IV-16) and (IV-17) are replaced by equations of the same form, but 
with the factor in in place of the operator d/dt 0 . Solution of these linear problems are available 
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[14, 2] and can be represented in the following form: 
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Here we consider that Re{\/z} > 0. 

The complex amplitudes of the heat fluxes can be determined as: 

Qmn = B n Tmn ^rnn> ^mn (^7 Pmn 'Lmn) ~b L) jnn 1 

H„ = l + /„ = ^LVScoth (A A) - 1 , 

Lti Ltv \ Ltv J 

Cmn = 2-~VinFm n (l), 
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(IV-23) 

(IV-24) 

(IV-25) 


(IV-26) 

(IV-27) 

(IV-28) 


3. Complex Amplitudes 

Substituting asymptotic series (IV-8) - (IV-15) and (IV-7) into (IV-1) - (IV- 4 ) and collecting terms 
for the same powers of e, we obtain the following linear inhomogeneous equations to determine 
the unknowns in the mth order approximation: 
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(IV- 32 ) 
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(IV- 33 ) 


where the right-hand sides, Y%\j = depend on the low-order approximations. 

To obtain equations for the complex amplitude of the nth mode of oscillations we can replace 
d/dt 0 in these equations with in. Consequently, equations (IV- 29 ) - (IV- 33 ) and (IV- 26 ) for the 
complex amplitudes can be represented in the following matrix form: 

(IV- 34 ) 
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(IV- 35 ) 


v£ = c mn , y£ = D m ». 

Let us now specify the structure of the right hand side terms Y mn . We can represent them as. 

Y — S 4- N + F mn (IV- 36 ) 

I mn — ^mn 1 ^mn ' 1 run 1 v 

where S mn is generated by slow time scale evolution of linear terms, N mn is generated by non- 
linear terms, and F mn is the external forcing. S mn can be found from linear equations (IV- 29 ) - 
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(IV-33) , (IV-16), (IV-17), (IV-24), (IV-25), (IV-27), and (IV-28) and asymptotic series (IV-7) and 
(IV-8) - (IV-15): 
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where the superscript in the brackets near S mn shows the number of the vector S mn component 
and 8ij is the Kronecker delta. 

In our case, F mn is non-zero only for n = 0,1. The order m for which F m o ^ 0 depends 
on the relation between the small parameters e and <5 determined by (III-33) and (III-37) . For 
S ~ e m * we have F m o ^ 0 only at m = m,. Since e and 6 are independent parameters, m, can be 
selected arbitrarily. However, a better way is to set m, using the characteristic asymptotic form. 
This can be done by the following physical reasoning, e and 8 are parameters responsible for two 
instabilities regarding a vapor bubble. The former is responsible for the instability due to rectified 
heat transfer and the latter for the instability due to the surface tension and deviation from the 
saturation state. The characteristic times required for development of these two instabilities are 
proportional, respectively, to the energy of the acoustic field, or e 2 , and to the deviation from the 
equilibrium state, or 6. These two instabilities can be brought to the same order of approximation 

if 


8 


~ e 


(IV-37) 


This is a condition for the characteristic asymptotic form. Accepting (IV-37) we have the following 
expressions for the non-zero components of F mn : 



(IV-38) 


An expression for N mn can also be derived using governing equations and asymptotic, repre- 
sentations. We do not reproduce them in general form here since they are too unwieldly. 
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4. Evolution in Slow Time Scales 


To obtain equations describing evolution of the unknowns in slow time scales, let us consider the 
operator Mq, 
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The determinant of 

this matrix 

is non-zero: 
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O 0 







From (IV-37) we have L c /a Q ~ <5 ~ e 2 < 1, or detM 0 = 0(e 2 ). According Cramer’s rule the 
zero-mode components of solution, are 


(i) _ det 
m0 det Mo ’ 


i = 


(IV-41) 


where is the matrix M 0 with the ith column replaced by the right hand side vector, Y m0 . 

By definition of an asymptotic expansion X^o should be of order of unity. Thus 

detM£, = 0(e 2 ), i = l,...,7. 

The form of matrix M 0 (IV-39) shows that this condition automatically holds for i = 2 ,..., 7. 
However for i = 1 it becomes non-trivial. Calculating det Mq X \ we have: 


r(7)\ 
mO J 


detM™ = -^M a (^-YS+y3)+ Sj ^-{yiS-YS->yi 

where Y^o and Y%o can be found as limits of (IV-24), (IV-25), and (IV-27) at n — *• 0 : 

00 /* 
y 2 = CmO = J vfmo{ri)dv, y 2 = D m o = J tfW*?)*?- 


(IV-42) 


To avoid secular terms in expansion (IV-8) , Omo and det M^J, should be limited at tj — ► oo, j 
1, 2, ... This is the only requirement for the uniform asymptotic expansion which is not unique. The 
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non-uniqueness of the asymptotic expansion can be easily demonstrated for the case of L a = 0 
(the effect of surface tension is neglected). In this case det Mo = 0, det Mq^ = 0, and a m o is 
determined as an arbitrary limited function of slow times. If we set cWnO = 0, for m = 1,2,..., 
then we uniquely select the expansion in which <2o is the bubble radius averaged over the period 
of oscillation, (a) , since (o) = do (1 d* £Oio + £ 2 &20 ...) . Such a definition of arbitrary a m o was 

used in a previous investigation [2]. Similarly, for the case L a / 0 we can specify 

(IV-43) 


detM^ = 0. 


According (IV-41) this leads to 

a m0 = 0, for m = 1,2, ..., (a) = a 0 , (TV-44) 

which is a convenient definition of Oq. 

5. First Order Approximation 

In the first-order approximation the non-linear term Nj n in (IV-36) is zero, the components of 
the external forcing are given by (IV-38), while the slow-scale evolution vector Si n has only two 
non-zero components (IV-37): 
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Using (IV-38) we can explicitly resolve (IV-34) for the first-mode using Cramer s rule. 
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j = 1, ...,7,where are the components of the vector X n and can -be represented as: 
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Note that these equations are consistent with previous results [30, 17, 3, 59, 2]. 

From (IV-42) and (IV-43) we have using (IV-36), (IV-38), (IV-45), and (IV-35): 

p=0, 0 „ = a„(( 2 ,i 3 ,...). (IV -47) 

ut\ 

To determine a slow-time evolution of the mean bubble radius we need to consider the second-order 
approximation. 


6. Equation for Rectified Heat Transfer in the Second 
Order Approximation 

Since ao does not depend on t u the same is true for all other variables. Thus, the non-zero 
components of vector S in the second order approximation are 
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(IV-48) 


The non-zero forcing term is (IV-38), while the components of the non-linear term N 2n can be 
found in the second-order approximation from quadratic non-linearities. The non-zero components 
include modes 0 and 2 which correspond to the average fluxes and nonlinear doubling of the 
oscillation frequency. To obtain the equation of rectified heat transfer it is sufficient to consider 
just zero-mode vector N 2 o- After some algebra one can find the components of this vector m the 
following form: 
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Here E(z ) is the third-order exponential integral, 
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and 7 r" is the dimensionless second derivative of the saturation pressure along the saturation curve. 
Using the Clapeyron-Clausius equation (III-19) for a perfect gas, we have 
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T 2 d 2 p s 
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(2 -p) 

k!(i -pf 


From (IV-36), IV-38), (IV-48), (IV-42), and (IV-49) we can obtain the following equation for 
rectified heat transfer: 
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W 2 (ao) = aoV(a 0 ) 


(a + ^) + AfeNff - ^ («g> - *«£>) + 
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This equation agrees in limiting cases with the equation for rectified heat transfer obtained using 
simplified models [37, 3, 17, 2j. Note that Equation (4.2) in Reference [2] contains misprints. 


7. Third Order Approximation 

It is noteworthy that the asymptotic scheme used does not provide uniformly valid expansions in 
orders higher than 2. The source of asymptotic singularity is connected with the infinite spatial 
region for the thermal problem in the liquid. Below we consider this problem in detail. 


a. Spatial Matching of Asymptotic Expansions 

The above procedure for solution of the thermal problem is correct for determination of the nth 
modes of fluxes when n ^ 0. At n = 0 the temporal derivatives with respect to to in the left-hand 
side parts of equations (IV-16) and (IV-17) are zero. In this case the temporal derivative is treated 
in a straightforward way which causes a singularity of the asymptotic expansion (small parameter 
near the highest derivative). In the slow time scales a temperature perturbation propagates from 
the bubble to the bulk of liquid to characteristic distances 7 ? 1, and the slow spatial scales 

should be considered. For the boundary value problem matching of inner and outer asymptotic 
expansions is an appropriate asymptotic procedure. 

The inner problem for u m o follows from (IV-16): 
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where tf m0 is the matching constant to be determined. Solution of this problem is 
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J (IV-51) 
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and for zero mode of the heat flux we have 
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Comparing with (IV-26) we can determine 

OO 

CmO = 'dmQ + J f}fmo{ T l)d r l- 


(IV-53) 


For simplicity we limit ourselves to the third order approximation requiring a one-term outer 
expansion. In this case it is enough to introduce one slow variable ry 2 = eq (otherwise we introduce 
for each time scale, t k , its own slow spatial variable, rj fc = e k/2 r], k = 2,3,... ) and consider 
a zero-order approximation, u ( 0 o) , to outer solution u (o) (otherwise we consider a series, u 0 
yjg) eu W 4 . ### ). The outer problem in slow scale (£ 2 ? 772 ) is therefore 
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Solution of this problem with the following initial and boundary conditions: 
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where $20 I s an arbitrary function of time. 

To match this solution with the inner solution, we substitute r? 2 = eq and expand this expression 
at € — *• 0 : 
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< 9 C y/^2 ~ C 

At the same time the two-term outer expansion of two terms of the inner solution, Uq , (IV-51) is. 


Uq ) = e 2 u 2 0 + e 3 U 3 o + — — c 2 < $20 + 


T 2 o — $20 + 


00 

J X (x - 1) / 2 o( 




+ £ 3 $30 + — 


(IV-57) 


Here we used the fact that / 30 = 0. This is true for harmonic oscillations, since the regular non- 
linearity of the equations generates modes 0 and 2 in the second order and modes 1 an me 
third order approximation. 
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Comparing (IV-56) and (IV-57) one can find: 


$20 = 0, 


$3° = — 7= 

yJ-K 


$20 = y~ (^20 + E 20 ) 

Lri 


t2 

$200 f 

yfo J 


$$20 


Vt 2 — C 


£20 = j x(x-l)f 20 (x)dx, 
1 

$200 = $20 |t 2 =0 • 


(IV-58) 


b. Equation for Rectified Heat Transfer 

In the third order approximation we have 


V W - 

I on — 


1 dap 


y( 2 ) _ n 

^30 “ u > 


v ( 3 ) __ 37 dap 

^30 “ ^ Ox > 

Oq OI 3 


r3 ° ” Oo&a’ 30 

ri 4) = y 3 ( 0 5) = 0, yjj> 6) = c 3 o = $30, yi 7 ) =o. 

General conditions (IV-43) and (IV-42) provide us with the following equation: 

^ = W 3 , W 3 = L 0 V{ao)ti 30 . 

% 

The zero mode of temperature in the second order approximation can be found using (IV-39)- 
(IV-41), (IV-36), (IV-48) and (IV-49): 

(IV-60) 


(IV-59) 


We also can determine 


T 2 o = C 20 + AZ? 2 o - N20 + 37 N io) ‘ 

a 0 Vi 


E 20 — “ ^- e 1 a nTu 


1 + 


Lri 




(IV-61) 


Equations (IV-58), (IV-60) and (IV-61) specify $ 20 as a function of ap. 

Generally, if we know equations in time scales t 2 , £3, ... then we can introduce one slow time, 
t — and represent the equation for rectified heat transfer in the form. 

da* = da 0 (£ 2 , t 3 , ■■■) ^qq _ = w / ( T) €T> ...)) + eW3 ( r , €r , oo (r, er, ...)) + - 

A- (IV-62) 


ch cft 2 ^3 

Thus, we have for ap the following equation: 


dap 

dr 


= W 2 (a 0 )- 


€LpV{(i 0) 


T 

$ 2 00 f d$2o( Q o) *9 Q 0 (C) $C 
/F J da 0 dC y/r-Q 


+ 0(6 2 ). 


(IV-63) 


The term W 3 reflects an important physical effect which shows that the growth rate depends 
not only on the current value of the bubble radius, but also on the history of the development of 
the thermal boundary layer. Mathematically the type of equation also changes from an ordinary 
differential equation in the second order approximation to a non-linear integro-differential equation 
in the third order approximation. 
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c. Thin Boundary Layers 

The case Ln ~ eao corresponds to a thin oscillating thermal boundary layer in the liquid and 
requires some special asympotic treatment. In this case the thickness of the ‘fast’ boundary 
sublayer is of order Ln/o, 0 <C 1. The evolution of the ‘slow’ thermal boundary layer in the liquid 
occurs at 77 — 1 Lri/ao, which, therefore, is not limited by the case 77 1, which we considered 

earlier to treat the outer problem. Instead for the outer problem we need to consider in the third 
order approximation a general case 77 ~ 1. 

Introducing the boundary layer variable, Q = do (77 — 1) /L77, we can represent the inner 
problem (IV-50) up to the third order approximation in the form 


d 2 u 20 

dQ 

U 2o|q. = o 

“30| C/=0 
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+ e 
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dQf dC f 


(o) 


C/=oo 


H 

d rj 


77=1 


where u ^ is the outer solution. Solving this problem and taking into account that /30 = 0 and 
T30 = 0 we can determine the zero-mode heat flux in the third order approximation (IV-52) 


Q30 


do chi30 
Lti 


= ~^30- 

</= 0 


(IV-64) 


Solution of these equations allows us to determine the following 
problem and the boundary conditions and the heat flux: 


relations between the outer 


u o | = ^20 + ■£'20, 

1,7=1 


^30 = 



+ 


dui o) 


drj 


T)=l, 


(IV-65) 


where the function £20 can be computed using (IV-61) since it is a limiting form of (IV-58) at 

Lri/a 0 1. 

Consider now the outer problem (the problem at 77 ~ 1). For simplicity we pass the mass 
flux term proportional to p 1 and the heat source due to viscous dissipation ( Ln/L d ) <C 1 to 
higher order approximations (see estimations (IV-19)). The outer problem will be then: 


c 2 Qq 

Lti 


du^ d\na 0 
dt 2 dt 2 




dr] 


1 d f 2 d u< o' 1 
rf drj \ ^ drj 


= 0, 


(IV-66) 


subject to the following boundary conditions: 


u, 


(o) 


,7=1 


= ^20 + E- 


20 , 


U , 
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(IV-67) 
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This problem can be solved with the following initial conditions: 


(o) 


t2=0 


= B(n), 


(IV-68) 


where H(t ] ) is some initial temperature profile near the bubble. For a bubble generated in a liquid 
of uniform temperature one can set H(rj) = 0. 

Generally solution of problem (IV-66)-(IV-68) can be obtained only numerically. This is the 
essence of a numerical-analytical method, which treats the problem in the third order approxima- 
tion and requires numerical solution of the boundary value problem (IV-66)- (IV-68). Numerical 
solution allows us to determine $30 in (IV-65) and evolve the bubble meanradius in a slow time 
scale according (IV-59) and (IV-62): 


dap 

dr 


= W% (ao) + eL /3 V (ao)l?30- 


(IV-69) 


Note that this method generalizes the consideration given above. If we drop the asymptotic 
assumption Ln/a 0 ~ e < 1 (consider Lri/a 0 ~ 1) then we obtain analytical solution (IV-63). 


8. Vapor- Gas Bubbles 

The basic asymptotic technique for vapor-gas bubbles is similar to that described above for pure 
vapor bubbles. There are several differences, which we address in this section. The major difference 
is in the heat and mass transfer in the gas where the energy and diffusion equations are strongly 
coupled. The energy and diffusion equations in the liquid are also coupled, but only through 
the boundary conditions. They also have similarity, which does not require special solution of the 
diffusion problem, and equations for the energy equation can be used (a small term representing the 
heat source due to viscous dissipation may be dropped). Therefore the complex amplitudes of the 
boundary temperature and concentrations can be determined by solving a general system of linear 
equations in the mth order approximation, while the zero-modes of the heat and mass fluxes in the 
liquid will be related to the boundary values by similar relations as for a pure vapor. Such system 
can be solved in the second and the third order approximations to obtain an evolution equation of 
type (IV-69). In the present study we, however, limit ourselves to the first-order approximation 
to evaluate the influence of the inert gas component on the bubble forced oscillations. 

a. Thermal and Diffusion Problems 

In the first order approximation we seek solutions of the governing equations m the form: 
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T = 

1 m 
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£rn = 
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e Re {q’ ma e 


**}, = m = l,g. 

1 0-0 


Liquid 

The heat and mass transfer in the liquid is described by equations 
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“ °- n '-' =r " ^ 

1 ( Lti\ 2 1 5 / 2 dc' t \ _ /i _ / 

ic] ; — ~ 1 V ~7T — U > C /lT 7 =l — c l 

1 Le t \ a 0 J V 2 dr] \ dq J 

Solution of these problems is well-known and can be found elsewhere [30]: 
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Gas 

The thermal and diffusion problems in the gas are coupled. It is convenient to represent them in 
the form 

r .u- r - (IV - 70) 


1C , - 


1 / Lr« 


Lcq V &o 


2 r 


LA \ 4- ^1— ( ri 2 ^—^\ 


rf dr ] \ dr] ) ?? 2 dq \ dr] J 


— 0) C iln=l C ia' 


Solution of these equations can be sought in the form 


1 , 1 , 
V - k~p' ~ - exp At;, q exp A 77 . 

9 V v 


The eigen values A can be found from the characteristic equation 


A 4 - i (Le g + 1 + k T k c ) 


J Tg 


do 


A 2 - Le a = 0. 


36 


(IV-71) 


This equation has four roots: 

A m = W = -■»-, rn = 1,2, 

V 2 ^Tg 

Sl>2 = i {(Le ff + 1 + ^r^c) ± [(M, - l) 2 + 2 k T k c [Le g + 1) + /c£fc c 2 ] 1/2 } . 
A non-singular solution in the domain 0 ^ tj ^ 1 can be written in the form. 


TL - kyp' = - [Ai sinh (Aip) + A 2 sinh (A 2 t?)] , 

9 ' T] 

d — — [Bi sinh (Aip) -I- H 2 sinh (A 2 tj)] , 

V 

where the constants of integration A m and B m can be found from equations and boundary condi- 
tions (IV-70). The heat and mass fluxes then can be found using relations (III-ll). The dimen- 
sionless complex amplitudes of these fluxes can be represented in the following form 


ia = CJc + (K - M) fr, C = c'o9c + (T' a - M) g T , 

. (1 - s 2 )n - (1 - si)r 2 
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Sl - S 2 


q m = [1 - s m + i {Le g + k c k T )\ Kn coth A m , r m = (1 - s m + ik c k T ) A m coth Xm, m- 1,2. 

These relations can be substantially simplified for small concentrations of one of the components, 
when we have kx — * 0 and, correspondingly . 


Sj L/€>g ) S 2 1? 


f c = l — coth Ai , 5t = 1- A 2 cothA 2 , It = 0, 

g c = [(1 - Leg) fc 4- iLe g (A 2 coth A 2 - Ai coth Aj)] . 

Leg (1 — Leg) 

Note that since g c ± 0 we have some influence of the concentration profile on the heat flux. 
However, this influence is small (since d a should be small at small c o0 ) and can be also neglected. 
So in this limiting case we can consider a decoupled system. 

b. Equations for Complex Amplitudes 

The set of linearized equations for complex amplitudes can be represented in the following dimen 
sionless form 

ia' - w' - p (£ 4- O = 
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^/31 
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^l(^) 4-c«£ + (i-<*>)£ = °' 
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(1 - p) L/3 


k v — 1 + p 

From a computational viewpoint, it is easier to solve a 8x8 linear algebraic system directly than 
to derive an analytical form. We tried two standard routines: L U decomposition and Gauss-Jordan 
[71]. Both work well and provide fast and accurate solution of the above system. 
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Chapter V. 

Asymptotic Technique For Standing 
Waves 

For description of the vapor bubble motion in a standing acoustic wave we use the same multi- 
scale technique as for bubbles in isotropic acoustic fields. We limit ourselves to the second order 
approximation. There are some peculiarities for standing waves which we consider below in detail. 


1 . Asymptotic Expansions 

In the case of standing waves we have additional unknowns compared to the case of vapor bubbles 
in isotropic fields, which we expand in the following asymptotic series: 


Xb {t) = (to, +...], (y-i) 

U b (t) = (V-2) 

m b {t) = ~np v <x>a-l [1 + (*o» h, •••) + •••] • ( V-3 ) 

We also have • 

sin {x 0 [1 + ex x + ...]} = sinx 0 cos [€XqX\ + ...] + cos xo sin [ex 0 xi 4- ...] 

= sin Xo 4- exoXi cos xo + ... 

cos {x 0 [1 + exi + ...] } = cos x 0 cos [ex 0 X! + ...] - sin x 0 sin [ex 0 Xj + ...] 

= cos Xo — exoXj sinxo + ... 

We can substitute these expansions together with expansions of variables used for description 
of bubble dynamics in isotropic fields into the governing equations and collect terms near the 
same powers of the small parameter £. This yields the following additional equations related to 
the bubble translational motion: 
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The matrix M n has a non-zero determinant, and solution, Xm n , can be obtained either analyt- 
ically, or (which is easier due to unwieldly analytical formulas) numerically by standard methods 
of matrix inversion. Note that the first seven equations form an independend subsystem (similar 
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to the system for bubbles in isotropic acoustic fields). If solutions Xml, j = 1,— .7 are found 
then the rest of unknowns can be determined as 


X (6) 

^ mn 

x v» 

x m 



l 

inx o 






Let us now be more specific about the structure of the right hand side vector. As before we 
can represent the right hand side terms Y mn in the form: 


Y mn = S mn 4- N mn + F mn , (V-6) 

where S mn is generated by slow time scale evolution of linear terms, N mn is generated bv^ non- 
linear terms, and F mn is the external forcing. Slow time scale evolution components, Smn, for 
k = 1, ..., 7 are the same as in (IV-37), while the other Sml are: 
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Uzing the same notation for F mn we can find that in our case only F m0 , Fn, and Fmn are 
non-zero. As we showed in the previos section to realize the characteristic asymptotic form we 
need to satisfy condition (IV-37). We also have the instability of the undisturbed state due to the 
gravity. To bring it to the second order approximation we accept the following condition: 


(V-7) 


2. Evolution in Slow Time Scales 

To obtain equations describing evolution of the unknowns in slow time scales let us consider 
operator Mq : 
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Note that the first 7 equations form a closed subsystem with a non-zero system matrix, Mo^s, 
having determinant 


det JVto,su(> — 


67 (1 - p) kjjp Lg 


Thus we can determine using Cramer’s rule: 
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All other equations including the slow-time scale evoluton condition for the 7x7 subsystem are 
similar to those obtained for isotropic fields (??)-(IV-43). 

The last two rows form of the system matrix Mo (see (IV-39)) correspond to equations 
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The third and the last condition for the right hand side vector can be obtained if we consider the 
eighth equation: 


-3u = *£. 


We have using (IV-41) and (IV-40): 
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det 'y 

Solutions of the system are not unique, and, as in the case with a m o, arbitrary x m o can be 
assigned. To specify a solution we set x m o = 0 for m = 1, 2, .... In this case we also have (x) = Xq. 


3. Linear Approximation 

In the linear approximation we have 
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We can explicitly resolve (III-39) for the first-mode using Cramer s rule. 
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where are the same as (IV-47) for isotropic fields. 

Then we can determine the complex amplitudes describing the bubble mass, position, and 


translational velocity: 
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4. Equation for Rectified Heat Transfer 


Since cio does not depend on t lf the same is true for all other variables. Thus, the non-zero 
components of vector S in the second order approximation are 
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The components of non-linear term N 2n (V-6) can be found in the second-order approximation 
from quadratic nonlinearities. The non-zero components include modes 0 and 2 which corresponds 
to the average fluxes and nonlinear doubling of the oscillation frequency. To obtain the equation 
of rectified heat transfer it is sufficient to consider just zero-mode vector N 20 . After some algebra 
one can find the components of this vector in the following form: 


n£? = SReK^n}, (V-17) 

N$ = 0. 

iV 20° ) = Re { " PCU ) UU } ‘ 

Components = 1, ...,7 are similar to those for isotropic fields (IV-49). 

Therefore from conditions of existence of periodical solutions (IV-42) and (V-12) we can obtain 
the following equation for rectified diffusion in a standing acoustic wave: 
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Chapter VI. 

Numerical Methods 


For computation of bubble dynamics in acoustic fields we applied the following three numeri- 
cal methods: a purely numerical integration of the detailed governing equations, the numerical- 
analytical method, and the analytical solutions described above. Each of these method has its 
own advantages and deficiencies, and the best method to use depends on the task performed. 
For example, for computation of bubble dynamics during relatively short times (of the order of 
hundreds of oscillation cycles) a purely numerical integration provides the best results in terms 
of accuracy and reasonable utilization of computational resources. This method, however, is not 
applicable for computation of bubble dynamics during the many thousands, and even millions 
of oscillation cycles required for bubble stabilization in an acoustic field. In this case the best 
results were obtained using the numerical-analytical method, which approximately utilizes the 
same computational resources as computations using analytical formulae. Despite the fact that 
the numerical-analytical method is much faster than the straight-forward numerical integration 
(computation of a test case of 3600 cycles of bubble oscillation required 30,000 times less CPU 
time for the numerical-analytical method than for the straight-forward method), it is still too slow 
for parametric studies such as the investigation of the equilibrium states, and the investigation of 
dependencies of bubble parameters in acoustic fields on multi-parameter inputs (say the accom- 
modation coefficient, liquid superheat, acoustic amplitude, and frequency). For this task the best 
method is the analytical method. 

We have developed several codes based on different models of vapor-gas bubbles. However here 
we describe the methods for pure vapor bubbles only which were the best tested and validated by 
comparisons with each other and with the methods used by other authors. 

1. Straight-Forward Finite Difference Scheme 

The system to be solved can be considered as a system of two algebraic equations and five time- 
derivative equations that include the mass conservation equation at the bubble surface (IV-1) to 
determine a, the modified Rayliegh-Plesset equation (IV-2) to determine w la , the vapor energy 
integral (IV-2) to determine p v , and the two energy equations, (IV-5) and (IV-6), for the vapor 
and liquid temperature, respectively. Among the five time-derivative equations, the first three are 
ordinary differential equations, while the last two are partial differential equations. To solve the 
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partial differential equations, one can apply a finite difference scheme for the spatial derivative 
terms in equations (IV-5) and (IV-6). The infinite computational domain required in equation 
(IV-6), however, may cause difficulty in numerical implementation. For numerical convenience, we 
can further transform equation (IV-6) by specifying £ = l/ry, such that the computational domain 
of equation (IV-6) is from 1 to 0. In terms of £ the liquid energy equation (IV-6) becomes 
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By applying the second order central differencing scheme for the spatial derivative terms, one 
can rewrite equations (IV-5) and (VI-1) as 
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If there are N v and Ni nodes used for discretization in equations (VI-2) and (VI-3) , respectively, 
then equations (VI-2) and (VI-3) can be considered as two sets of ordinary differential equations 
in the u nk nowns. Arranging a,w ia ,p v , and all T V i and Tu into a single vector of- unknowns X, we 
can write equations (IV-l)-(rV-2), (VI-2), and (VI-3) in the form: 


X = F(X,t). (VI-4) 

The vector function F is shorthand for the right hand sides of the equations. A fourth-order 
Runge-Kutta scheme is applied to integrate equation (VI-4) over time. 

To start the time integration, the initial condition for the temperature distribution inside the 
bubble and in the liquid need to be specified. Usually we assume that the vapor and liquid 
temperature are initially uniform, 


HU = 21U = r». (vi- 5 ). 

During the time integration the boundary conditions for the temperature are satisfied by enforcing 
the zero flux at 77 = 0, and the temperatures at the interface and far from the bubble: 
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(VI-6) 


r.U = n| ( „ = T„, T, | (rf = 

To satisfy the zero flux boundary condition we specify T v \ = T v2 . 

To obtain the bubble wall temperature, T a , we solve two algebraic equations (IV-3) and (IV-4), 
where we use the following approximations for the saturation pressure and temperature consistent 
with the Clausius-Clapeyron equation: 
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Since (IV-3) and (IV-4) are nonlinear equations, an iteration procedure to determine T a is required 
using p v , q va , and qi a obtained from the previous time step. We use the following converging scheme: 


° E 2 {T£)' 


(VI-7) 


where the superscript shows the number of iteration (for initial approximation we use the value 
of T a from the previous time step), and functions Ei(T a ) and E 2 (T a ) are calculated using the 
following formulae: 
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Since the thermal boundary layers on both sides of the bubble wall are very small and the 
gradient of temperature is usually large, a very fine grid must be created near the bubble wall on 
both sides to adequately resolve the thermal boundary layers. To avoid the use of large number 
of grid points in the computational domain, an unequal-spaced grid which has highly clustered 
points near the bubble wall on both sides is applied. Figure VI-1 shows a typical temperature 
distribution during the bubble oscillation and the grid distribution along the space. 
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Figure VI-1: Typical temperature profile inside (77 < 1) and outside (77 > 1) a vapor bubble 
oscillating in an acoustic field. The dots on the curve correspond to nodes of the computational 
grid. 


2 . Numerical- Analytical Method 

In this section we describe a numerical algorithm suitable for solution of the linear partial differ- 
ential equation (IV-66). For convenience, we rewrite this equation in the form, 


where 
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The conversion, u = qu^ transforms this equation into the one-dimensional heat conduction 
equation, 


dt \&n vj 

with the same boundary conditions given by equation (IV-68). 

We discretize time and the radial coordinate and approximate the first and the second spatial 

derivatives with second order of accuracy: 
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(VI-10) 
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Applying the first order approximation to the time derivative and writing equation (VI-9) in 
the implicit manner results in a system of algebraic equations of the form, 

where 


where: 
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o' = Si kt/s + Gip m At, 

7 = S^Ai/s — G^PmAt, 

(3 = 1 + a + 7 + p m &t/r\ m . 

The tridiagonal system (VI-12) can be solved effectively by using a standard algorithm [71]. 
It assumes that the values of function u at the neighbouring nodes are linearly dependent, 


w m_l — Am-l u m + B m - J. 


(VI-13) 


Substituting equation (VI-13) into (VI- 12 ) one can find that the coefficients An, B m satisfy 
the following recurrent relations: 


_ a D = < 1 +7-Sm-l (VI-14) 

Am ~ P-^Am-l m P-lAn-l ' 

We assume that the temperature at the bubble boundary is presribed, i.e., A\ = 0, B\ = U a {t n ), 
while at “infinity” the heat flux is zero, i.e. u M = B M f (1 - Am)- Then, by applying successively 
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equations (VI- 14) one can find the coefficients A , B and then back substitute them to find the 
sought function u. The outlined numerical procedure is absolutely stable because the coefficients 
A m are less then unity, {Ami < 1. This statement can be proven by deduction, since A\ = 0 and 
P — 7 = 1 + a > a. 

The index M corresponds to the right boundary of the computational domain. The size of the 
computational domain expanded in time as the thickness on the bubble’s boundary layer grows 
up. We used the following algorithm to track this expansion. First, we specified the dimension size 
L of the array coresponding to temperature, u = T/rj, and we set, M = L/ 2, initially. Second, we 
computed the temperature distribution using equations (VI-13) and (VI-14) for the index range: 
1 < m < K, where K = M + 10. Then, we computed the spatial derivative duj dr) and found 
index J such that M < J < K and \{du/dr))j\ < 1CT 8 . We updated the value of index M to J 
on the new time step. This procedure was repeated until index M reached the value of L — 10. 
When this happened the grid size was doubled and a cubic spline fit was applied to compute the 
temperature distribution corresponding to the new grid. 

The dynamics of the mean bubble radius was then computed using the fourth-order Adams’ 
extrapolation scheme. When solution for the average bubble radius was found, we computed the 
current bubble radius by superposition of the third-order solution for the mean bubble radius and 
the linear solution for the fast oscillating part: 


a[t ) = Oo(t) (1 +eRe(a I1 e“ Jt )). 


(VI-15) 


3. Analytical Solutions 

The analytical solutions were used mostly for investigation of the parameter space. Computation 
of the complex amplitudes of oscillation were performed analytically for the case of a pure vapor 
bubble and numerically, using LU-decomposition, for vapor-gas bubbles. The results for vapor-gas 
bubbles were verified by comparing the limiting case of pure vapor bubble with the case of vapor 
bubbles and with physical interpretation of the computed patterns. The results for vapor bubbles 
including the first-order and the second-order approximations were validated against limiting cases 
obtained earlier [2] when the surface tension, liquid viscosity and compressibility were neglected. 
To plot the phase portraits of the autonomous systems of type (V-18) were used and integrated 
with the Runge-Kutta 4th-order method. 

For computation of the mean radius evolution in the third-order approximation nonlinear 
integro-differential equation (IV-63) was integrated using a Runge-Kutta 4th-order scheme with 
computation of the right hand side using an iterative scheme and evaluation of integrals based on 
a Gauss-type quadrature for modified Chebychev polynomials C n {x) (polynomials orthogonal on 
[0,1], not [-1,1]): 
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The weights and nodes of the Chebychev polynomials were found using standard procedures of 
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orthogonalization (e.g., [71]). In most cases 10 weights and nodes provided sufficient accuracy 
(tested by comparison with exact solutions for relative error below 10 6 ). 

After determination of the mean radius evolution, the current bubble radius was found by 
superposition of the first-order linear oscillation and the mean bubble radius according (VI-15). 

4. Comparisons of Bubble Dynamics Computed by Vari- 
ous Methods 

For validation of the numerical methods, we compared them with each other and with published 
computational results on vapor bubble dynamics in acoustic fields from Hao & Prosperetti (1999) 
[35] , who used a straight-forward pseudospectral numerical method. The model used in Ref. [35] 
ignores the nonequilibrium phase transitions and is applicable for a low frequency field. To obtain 
a ‘quasiequilibrium’ scheme in our model, we set /3 = 10 6 (for the computed case of 1 kHz acoustic 
field and millimiter size bubbles the difference between computations with (3 = 1 and (3 = 10 
were very small). Comparisons between all methods are shown in Figures VI-2 and VI-3. 

The straight-forward numerical computations fall on top of the cited results [35]. The other 
two methods give results within 10% of relative error for maximum deviation with the straight- 
forward computations. Note that much smaller error (about 2%) was obtained for the mean bubble 
radius computed using the numerical-analytical method and the straight-forward methods. This 
is related to the use of formula (VI-15) which neglects the second and the higher harmonics of 
oscillations, and may be improved by inclusion of the second and the third-order terms in (VI-15). 
The lower curve presented in Fig. VI-2 was obtained using initially nonuniform (quasisteady, or 
hyperbolic) temperature profile in the liquid, which demonstrates the importance of selection of 
proper initial conditions. 

Figures VI-4 and VI-5 show that agreement between the results obtained using the numerical- 
analytical method, the analytical solution and the straight-forward method is good. 
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Figure VI-2: Comparison between computations based on the asymptotic theory with the initial 
quasisteady temperature profile in the liquid (the lower solid curve), initial uniform temperature 
profile in the liquid (the upper solid curve), and the numerical results of Hao & Prosperetti 
(1999) (the dashed curve), which are on top of our results using the purely numerical method and 
quasi-equlibrium scheme of phase transitions. The bubble radius is normalized with the primary 
resonance radius, 2.71 mm. The initial mean radius in the computations using the present theory 
is 0.1 mm. 



Figure VI-3: Comparisons of computations using the purely numerical method (the dashed curve) r 
numerical-analytical method (the thick solid curve), and the third-order analytical solution (the 

thin solid curve). 
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Figure VI-4: Comparison between computations using the third order asymptotic theory with the 
initial temperature jump ($200 = $20 ( a in ) , a in — 98 pm; the solid curves) and purely numerical 
simulations with the initial bubble radius 100 pm using the detailed equations (the gray region 
and the dashed curve). Letters L, M, and U near the curves relates to the lower, mean, and upper 
slowly changing bubble radius. The curves L and U were computed by addition and substraction 
of the amplitude of bubble oscillation predicted by the linear theory. The numerical-analytical 
method gives results close to the straight-forward computations. 



Figure VI-5: Comparison between the purely (the dashed curve) and the numerical- analytical 
method (the solid curve). The other notations and computation parameters correspond to the 
case shown in the previous figure. 
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Chapter VII 


Analysis of Linear Bubble Dynamics 
in Acoustic Fields 

1. Vapor Bubbles 

Computations for water and helium of vapor bubble dynamics in acoustic fields were carried 
out for a range of bubble sizes, frequencies, and amplitudes appropriate for the present theory. 
Since water and helium vapors deviate from perfect gas behavior (particularly true for helium 
at low temperatures) the property values given in Table 1 were utilized in computations. Other 
quantities, such as the gas constant and the specific heat ratio, were derived using the perfect gas 
relations. 


Table 1. Properties of water and helium used in comp utations 


Parameter 

Unit 

Water 

Helium 

Too 

K 

373 

4.2 

Ps(Too) 

kPa 

101 

100 

Pvoo 

kg/m 3 

0.597 

16.3 

pi 

kg/m 3 

958 

125 

cpv 

kJ/(kg-K) 

2.03 

2.2 

Cl 

kJ/(kg-K) 

4.22 

6.76 

A, 

mW/(m-K) 

24.8 

2.3 

A i 

mW/(m-K) 

680 

27.1 

pi 

// N-s/m 2 

279 

3.57 

c 

km/s 

1.54 

0.18 

l 

km 2 /s 2 

2.26 

0.0209 

<r 

mN/ra 

58.9 

0.1 


Figures VII-1 and VII-2 demonstrate typical dependences of the amplitude, |an| , and phase, 
arg(an), of bubble radius oscillations on the average bubble radius, Oq. Computations were per- 
formed using (IV-46) and (IV-47) for water at 1 atm and an acoustic frequency of 60 kHz. Qual- 
itatively the dependence is the same for water at different pressures and frequencies or for other 
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liquids. Two limiting cases exist: the case of equilibrium (or more precisely quasi-equilibrium) 
phase transitions where we assume that the vapor pressure is equal to the saturation pressure at 
the interface temperature, and the case of the absence of phase transitions where ft = 0. The equi- 
librium case can be formally obtained if we set (3 = oo despite the actual value of P cannot exceed 
1. The computed curves show that the linear response of the vapor bubble to acoustic excitation 
strongly depends on (3. The two limiting cases determine two resonances. For (3 = 0 we have only 
the primary resonance due to the liquid inertia and the vapor elasticity typical for gas bubbles. 
In the case of equilibrium phase transitions an additional low-frequency ( “second” ) resonance due 
to phase transition and surface tension takes place. For water vapor bubbles at moderate {3 the 
amplitude of oscillation at the primary resonance is smaller than at the second resonance. For 
small P the bubble response curves show strong primary resonance. Note that at p = 0 formally 
there exists the “second” resonance at a w 2cr/(3p s ). However to reach this resonance we need to 
have A p = —3 p 3 which corresponds to negative liquid pressure. This fact was discussed earlier by 
Khabeev [30]. Moreover, results for small bubble sizes below the theory limit line in Figures VII-1 
and VII-2 violate assumption (III-37) and are not justified by the present theory. 

To obtain more insight into the nature of the vapor bubble resonances, it is convenient to 
represent the bubble oscillator in the form: 


a! + pa ' + ula' = -p', (VII-1) 

where p is the dissipation coefficient, u> n is the bubble natural frequency, and a' and p' are the 
nondimensional bubble radius and driving pressure, p and u n depend on the mean bubble size, 
do, the driving frequency, the accommodation coefficient and other properties of the liquid and its 
vapor. The natural frequency oj n can be found from equation (VII-1) by substituting there the 
expression for the driving pressure and the linear solution for the bubble radius. Fig. VII-3 shows 
uj n (ao; P) • It is seen that the natural frequency of vapor bubbles is very different from that of gas 
bubbles, and it is strongly influenced by p. At very small P a vapor bubble behaves similar to a gas 
bubble, while at relatively large P (in the computed case P > 0.1) the bubble natural frequency is 
below the driving frequency for any bubble size, and there are no resonant sizes at all (for smaller 
acoustic frequencies there are two resonant sizes for large P also). For bubbles smaller than some 
critical size, we have uj„ < 0 which is an effect of surface tension specific for bubbles with phase 
transitions. Such bubbles are exponentially unstable in the absence of external forcing, while they 
can be stabilized by an acoustic field and behave similar to the bubbles of resonant frequencies 
much smaller than the driving frequency. 

Figure VII-4 illustrates the temperature profiles inside and outside the bubble at a fixed mo- 
ment of time (to = 27rn). Computations were made for 50 pm bubbles using analytical solutions 
(IV-23) and (IV-23). All parameters were the same for all plotted curves, except for the value of 
the accommodation coefficient. It is seen that in the illustrated cases, P substantially influences 
the temperature distribution in the vapor. The influence of p on the temperature in the liquid 
is less, and there is no visible difference between computations with P = 0.04 and P = oo. At 
smaller /?, the temperature gradients in the vapor are much higher than those predicted by the 
quasi-equilibrium theory, and the heat flux in the vapor can be comparable with the heat flux in 
the liquid. This affects the interface temperature and the temperature profiles in the liquid. 
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Figure VII-1: Relative amplitude of the forced vapor bubble radial oscillation in a 60 kHz acoustic 
field. The numbers near the curves show the values of the accommodation coefficient, ( 5 . The curve 
marked as “equilibrium” is computed using the quasi-equilibrium scheme of phase transition. 



Figure VII-2: Phase shift between the driving pressure and forced radial bubble oscillation. No- 
tations are the same as in the previous figure. 
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Figure VII-3: The vapor bubble natural frequency vs the bubble mean radius in an acoustic field at 
various values of the accommodation coefficient. The curve c e’ correponds to the quasi-equilibrium 
scheme of phase transitions. 



Figure VII-4: Temperature profiles predicted by the linear theory at various values of the accom- 
modation coefficient shown near the curves and the quasi-equilibrium phase transition, marked as 
“equilibrium” . Profiles are computed for a vapor bubble of radius 50/zm at the phase of oscillation 
corresponding to t == 0.. 
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2. Vapor-Gas Bubbles 


Computations of the complex amplitudes of forced bubble oscillation were performed for a water- 
air system at atmospheric pressure and temperatures close to the water boiling point. Parameters 
varied include the frequency of the acoustic field, bubble size, concentration of the inert g as inside 
the bubble, concentration of the dissolved gas, vapor and inert gas accommodation coefficients, 
and the thermal diffusion coefficient. The basic properties of the system are given in the Table 2 
below. All other properties were derived using relations for perfect gas mixtures. The values of 
C pg at different concentrations were found by linear interpolation between limiting values for pure 
components. The values of K g and X g were assumed constant, since for the given system these 
quantities experience relatively small variations at changing concentration. 


(The 


Table 2. Properties of water-air system used in computations _ 


ie empty cens 
Parameter 

3 CUIItJbpUIlU IA 

Unit 

J biic vent 

Liquid 

Vapor 

Air 

Vapor-Air 

Interface 

Too 

K 

373 

373 

373 

373 

373 

Pa (Too) 

kPa 





101 

Dl, Dg 

mm 2 /s 

0.011 



37.8 


K h Kg 

mm 2 /s 

0.168 



20.4 


A/, \ 

mW/(m-K) 

680 



24.8 


Rv,Ri 

kJ/(kg-K) 


0.455 

0.29 




kJ/(kg-K) 


2.03 

1.01 



Pi 

kg/m 3 

958 





P'1 

fi N-s/m 2 

279 





c 

km/s 

1.54 





lv 

km 2 /s 2 





2.26 

li 

km 2 /s 2 





0 

H 

GPa 





6.73 

dH/dT 

GPa/K 





0 

a 

mN/m 





58.9 

dcr/dT 

mN /(m • K) 





0 


Figures VII-5 and VII-6 demonstrate the frequency dependence of the normalized amplitude 
and phase of the bubble radius oscillation. Various curves in the figure correspond to various 
concentrations of the inert gas inside the bubble with other fixed parameters (the coefficients 
k c and k T vary since they depend on the content of the gaseous phase). The case of c«> - 1 
corresponds to the case of a pure gas bubble, and Cj 0 = 0 to the case of a pure vapor bubble. 
The latter limiting case was verified for consistency with the previous computational results. At 
Cjo = 1 strong primary resonance is observed. The amplitude of the bubble oscillation at the 
primary resonance frequency decreases at decreasing concentration of the inert gas. In the range - 
of inert gas concentrations between 1 and 10% an interesting effect of low-frequency resonance 
(or, rather, instability) is observed. As can be seen in the figure, the amplitude of a 10 fan bubble 
oscillation drastically increases at do = 0.04 and frequencies of order 1 kHz. This resonance is 
controlled by the surface tension and the two-component nature of the bubble content. 
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This resonance was first reported by Nagiev & Khabeev [65]. They explained it by the qualita- 
tive difference between the low-frequency response of the vapor bubble, which oscillates in phase 
with the forcing pressure (due to the effects of surface tension and evaporation-condensation main- 
taining the vapor pressure on the saturation line), and the bubble of non-condensable gas, which 
oscillates in anti-phase to the forcing pressure (due to stiffness of constant mass gas). The bubbles 
containing a mixture of non-condensable gas and vapor can oscillate in any phase in between 0 
and 7 r with the external low-frequency field depending on the concentration. At some critical con- 
centration [65] (in our case Cio « 0.04) a switch occurs between these two qualitatively different 
cases which causes the type of resonance behavior shown in the figure. This mechanism is different 
from the condensation-evaporation ‘second’ resonance for pure vapor bubbles [32, 34, 35]. 

Note that at very low frequencies the gas-vapor bubble will oscillate in phase with the forcing 
pressure due to dissolution of the inert gas. Theories which neglect gas diffusion in the liquid 
[65, 66] provide different low-frequency limits for the phase and amplitude of the bubble oscillation. 
We compared our computations with published results of other authors [65] and found qualitative 
agreement for the range of frequencies over which diffusion in the liquid is negligible. We also 
found some quantitative differences, which can be explained by the differences in the models used. 

Figures VII-7 - VII-10 show the dependencies of the amplitude and phase of the bubble radius 
oscillation on the mean bubble size for a fixed frequency. Computations were performed for 
different contents of the gaseous phase. In sill cases the ‘second’ resonance is observed in addition 
to a strong primary resonance. Genersilly, the second resonance is caused by the effects of surface 
tension and phase transitions. In the limiting cases of Cjo = 1 and Qo = 0, phase transitions are 
due to dissolution and condensation-evaporation of one-component gas (or vapor). The phase 
transitions of a two-component mixture are more complex and lead to anomalous behavior of the 
‘second’ resonance as discussed above. The characteristic times of condensation are controlled by 
the liquid thermal diffusivity and times of dissolution are controlled by the liquid mass diffusivity 
which are usually of different orders of magnitude (Lewis numbers of liquid are usually large, Le t 3> 
1). This creates two distinct size/frequency scales: 1) a<) < Le t Lnri and 2) Oo 2> Le l Lri, 
where the inert gas can be considered as soluble or non-soluble, respectively. 

Effects related to gas dissolution can be important in the low-frequency/small-radius range 
a 0 < LeJ^Lji. In range ao Le\ X ^Lxi we do not expect such parameters as the inert gas 
accommodation coefficient, the Henry constant, and the concentration of the dissolved gas to 
influence bubble oscillations. To demonstrate this point we performed computations at different 
Pi and H with other parameters fixed for different frequencies (see Figures VII-11 and VII-12). 
First, we found that the effect of non-equilibrium dissolution of the gas can be observed only 
at extremely low Pi . Already Pi = 10~ 6 produces results close to the quasi-equilibrium phase 
transition of the inert gas (this case formally correspond to Pi — oo) . However, the models of 
non-soluble gas and soluble gas show a difference in the low frequency range ao Lei'Ln- 
If the values of Pi below 10~ 6 (which means that less than one-millionth part of gas molecules 
hitting the interface experience phase transition) are not of interest, within the frequency range 
Oq < Le[ X ^Lj'i (and therefore everywhere) the quasi-equilibrium scheme of phase transition for the 
inert gas is acceptable. The effect of the Henry constant on the bubble oscillations is ^appreciable 
in frequency range ao < Le~[ 1/2 Lri and negligible for frequencies of range a 0 > Le t Lti . Note 

that for 10 /im bubbles the effects of gas dissolution are negligible at frequencies of the order 100 
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Hz. The same is true for larger bubbles and higher frequencies. Therefore a simplified model 
neglecting the effects of gas diffusion in the liquid and dissolution is applicable for a wide range 
of bubble sizes and frequencies. For slow processes, such as rectified diffusion, the diffusion and 
phase transitions of the inert gas are important in slow time scales, but can be neglected in the 
fast time scale that is inversely related to the frequency of the field. 

While the effect of nonequilibrium phase transition for the inert gas is negligible, the effect 
of the nonequilibrium phase transitions for the vapor component is strong in high frequency 
acoustic fields. Similarly, for a gas-vapor bubble with appreciable vapor component, the vapor 
accommodation coefficient influences the bubble dynamics (Figure VII-13-VII-14). Although the 
value of the accommodation coefficient cannot exceed 1, we made computations using very high 
values of this coefficient to determine if the scheme of quasi-equilibrium phase transitions (the 
bubble interface temperature is equal to the saturation temperature) is applicable. It is seen that 
such computations indicated as ‘equilibrium’ are close to the computations with 0 V = 1. Curves for 
/3 V = 0.1 substantially differ from the equilibrium curves. Therefore the quasi-equilibrium scheme 
of phase transitions is applicable if /?„ is about unity. For very small radii, the effects of the inert 
gas dissolution create a resonance discussed above. Such influence of the vapor accommodation 
coefficient on the character of forced vapor-gas bubble oscillations suggests that measurements of 
this coefficient can be performed not only with pure vapor bubbles, but with vapor-gas bubbles 
with substantial inert gas content. 

It is generally accepted that the nonequilibrium character of condensation/evaporation should 
be taken into account only for high frequency fields (e.g.[35]). This also depends on the bubble 
size, and, as we found, on the concentration of the inert component. As an example we show in 
Figures VTI-15 and VII-16 computations of a 10 /in bubble response to acoustic perturbations over 
a wide range of frequencies with concentrations of the inert gas near the critical value (c*o = 0.05). 
It is seen that the second resonance frequency and the bubble complex amplitude at frequencies of 
about 1 kHz depend substantially on the value of the vapor accomodation coefficient for (3 V < 0.1. 

The dependence of the amplitude and phase of vapor-gas bubble oscillation on the accommoda- 
tion coefficient can be used for determination of the value of f3 v . For efficient diagnostic procedures 
the bubble size, the frequency of the acoustic field, and the concentration of the inert component 
should be selected from the range corresponding to the highest sensitivity of bubble dynamics to 
the value of this coefficient. An example shown in Figures VII-17 and VII-18 demonstrates depen- 
dencies of the amplitude and phase of forced bubble oscillations on 0 V at different concentrations 
of the dissolved gas, which can be selected arbitrary for a nonequilibrium vapor-gas bubble. It is 
seen that for /?„ < 0.1 the parameters of oscillations strongly depend on (3 V for small and moderate 
0 , 0 - Larger /3 V can be determined from the bubble dynamics at relatively small c*o since at higher 
Cjo the parameters of oscillations are not influenced by the value of 0 V . It is interesting that the 
maximum slope of the curves for the amplitude of oscillation at 0 V ~ 1 is realized at intermediate 
values of c#) (in our case c,o ~ 0.1), which indicates that addition of 10% of inert gas can be 
beneficial for diagnostics of (3 V . 

If @ v < 0.1, then its detection can be performed using the dynamics of bubbles of various 
sizes. However, for detection of larger values of 0 V smaller sizes provide better resolution. The 
curves shown in Figure VII-19 and ?? demonstrate that in the specified range of parameters, the 
detection of (3 V for water can be problematic for bubbles of radius 30 (im and higher if 0 V > 0.4. 
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Figures VTI-21 and VII-22 show that for the best resolution in determination of relatively large 
/ 3 V , the frequency of the acoustic field should be properly selected. Both high frequency and low 
frequency curves show a weak dependence on (3 V in region /3 V > 0.1, while intermediate frequencies 
of the order 50 kHz provide the highest sensitivity to the value of the accommodation coefficient. 

The influence of the Soret-Dufour, or thermodiffusion, effect on bubble forced oscillations in 
an acoustic field can be determined by variation of the thermal diffusion ratio, kx, or related 
thermodiffusion coefficient k c (III-12). Figures VII-23 and VII-24 display dependencies of the 
amplitude and phase of bubble radius oscillation on the mean bubble radius computed at different 
values of k c . From (III-13) we have k c = 1.37. Variation of k c over a moderate range (limiting case 
k c = 0 corresponds to the absence of the Soret-Dufour effect) shows its appreciable influence on 
the vapor-gas bubble response, especially for sizes near the primary resonance value. One possible 
applications of this observation is measurement of the thermodiffusion coefficient using vapor-gas 
bubble dynamics. 

To demonstrate the Soret-Dufour effect the dimensionless temperature and concentration pro- 
files realized inside an oscillating bubble at times ut = 27rn are shown in Figures VTI-25 and 
VII-26. The bubble size is selected to be close to the primary resonance at the given frequency. 
The computations were performed for different values of k c . The results show that larger gradients 
of temperature and concentration inside the bubble are realized due to the Soret-Dufour effect. 
The amplitude of concentration oscillation near the bubble surface increases with k c . This may 
explain the increase of the amplitude of bubble oscillation because larger amounts of inert gas 
near the surface make phase transition more difficult and increase the stiffness of the bubble. 
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Figure VII-5: Relative amplitude of the forced vapor-gas bubble radial oscillation for 10 [i m bubble 
in water as a function of the acoustic frequency. The numbers near the curves show the value of 
the inert gas concentration inside the bubble. 



log io(to/2n) 


Figure VII-6: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation for 10 fx m bubble in water as a function of the acoustic frequency. The numbers near 
the curves show the value of the inert gas concentration inside the bubble. 
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Figure VTI-7: Relative amplitude of the forced vapor-gas bubble radial oscillation in 1 kHz acoustic 
field as a function of the bubble radius. The numbers near the curves show the value of the inert 
gas concentration inside the bubble. 



Figure VII-8: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation in 1 kHz acoustic field as a function of the bubble radius. The numbers near the curves 
show the value of the inert gas concentration inside the bubble. 
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Figure VII-9: Relative amplitude of the forced vapor-gas bubble radial oscillation m 60 kHz 
acoustic field as a function of the bubble radius. The numbers near the curves show the value of 
the inert gas concentration inside the bubble. 



Figure VII-10: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation in 60 kHz acoustic field as a function of the bubble radius. The numbers near the curves 
show the value of the inert gas concentration inside the bubble. 
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Figure VII-11: Relative amplitude of the forced vapor-gas bubble radial oscillation for 10 /zm 
bubble in water as a function of the acoustic frequency. Computations are performed for different 
values of the inert gas accommodation coefficient, /%, and the Henry constant, H (the values are 
shown near the curves). The thick curves marked ‘equilibrium’ correspond to quasi-equilibrium 
phase transition of the inert component. 
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Figure VI 1-12: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation for 10 jivn bubble in water as a function of the acoustic frequency. Computations are 
performed for different values of the inert gas accommodation coefficient, and the Henry con- 
stant, H (the values are shown near the curves). The thick curves marked ‘equilibrium’ correspond 
to quasi-equilibrium phase transition of the inert component. 
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Figure VII-13: Relative amplitude of the forced vapor-gas bubble radial oscillation in 60 kHz 
acoustic field as a function of the bubble radius. The numbers near the curves show the value of 
the vapor accommodation coefficient, (3 V . The curves marked ‘equilibrium correspond to quasi- 
equilibrium phase transition of the vapor. 



Figure VII-14: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation in 60 kHz acoustic field as a function of the bubble radius. The numbers near the 
curves show the value of the vapor accommodation coefficient, j3 v . The curves marked ‘equilibrium 
correspond to quasi-equilibrium phase transition of the vapor. 
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Figure VII-15: Relative amplitude of the forced vapor-gas bubble radial oscillation foi 10 /xm 
bubble in water as a function of the bubble radius. The numbers near the curves show the value 
of the vapor accommodation coefficient, The curves marked ‘equilibrium correspond to quasi- 
equilibrium phase transition of the vapor. 



Figure VII-16: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation for 10 /xm bubble in water as functions of the bubble radius. The numbers neai the 
curves show the value of the vapor accommodation coefficient, 0 V . The curves marked ‘equilibrium 
correspond to quasi-equilibrium phase transition of the vapor. 
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Figure VII-17: Relative amplitude of the forced vapor-gas bubble radial oscillation for 20 /un 
bubble in 60 kHz acoustic field as a function of the vapor accommodation coefficient. The numbers 
near the curves show the value of the inert gas concentration. 



Figure VII-18: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation for 20 fMm bubble in 60 kHz acoustic field as functions of the vapor accommodation 
coefficient. The numbers near the curves show the value of the inert gas concentration. 
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Figure VII-19: Relative amplitude of the forced vapor-gas bubble radial oscillation in 60 kHz 
acoustic field as functions of the vapor accommodation coefficient. The numbers near the curves 
show the mean bubble radius. 



Figure VII-20: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation in 60 kHz acoustic field as functions of the vapor accommodation coefficient. The 
numbers near the curves show the mean bubble radius. 
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Figure VII-21: Relative amplitude of the forced vapor-gas bubble radial oscillation for 20 fxm 
bubble as a function of the vapor accommodation coefficient. The numbers near the curves show 
the frequency of the acoustic field. 



Figure VII-22: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation for 20 /xm bubble as a function of the vapor accommodation coefficient. The numbeis 
near the curves show the frequency of the acoustic field. 
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Figure VTI-23: Relative amplitude of the forced vapor-gas bubble radial oscillation in 60 kHz 
acoustic field as a function of the bubble radius. The numbers near the curves show the value of 
the thermodiffusion coefficient k c . 





Figure VII-24: Phase shift between the driving pressure and the forced vapor-gas bubble radial 
oscillation in 60 kHz acoustic field as a function of the bubble radius. The numbers near the curves 
show the value of the thermodiffusion coefficient k c . 
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Figure VII-25: Dimensionless perturbations of temperature profiles inside 50 Jim vapor-gas bubble 
oscillating in 60 kHz acoustic field. Profiles correspond to times ut = 27m. The numbers near the 
curves show the value of the thermodiffusion coefficient k c . 



Figure VII-26: Dimensionless perturbations of inert gas concentration profiles inside 50 /rm vapor- 
gas bubble oscillating in 60 kHz acoustic field. Profiles correspond to times ut = 27rn. The numbers 
near the curves show the value of the thermodiffusion coefficient k c . 
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Chapter VIII. 


Analysis of Nonlinear Bubble 
Dynamics in Isotropic Acoustic Fields 

1. The Equilibrium Radius and Its Stability 

There exists a possibility of stabilization of the vapor bubble radius in an acoustic field. Complete 
stability analysis cannot be performed based on asymptotic theory, since the residual terms, such 
as in equation (IV-63) can cause a growth or decay of bubble radius in slower time scales than 
that taken into account (secular terms). However, if the steady radius exists, then it is stable in 
any slow time scale. If it does not exist then there is no stable radius in low order approximations. 
That is why low order approximations can provide valuable information on the bubble stability. 
Everywhere below when we speak about stability we mean stability within the framework of the 
third-order theory, which can potentially be violated in slower time scales. 

Let us define the steady radius, a*, as a zero of function W 2 : 

W 2 (a,) = 0. (VIII-1) 

Consider a small perturbation of ao near this radius: 

a = ao -a.. (VIII-2) 


Linearizing (IV-63) near ao = a, we find: 
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Solution of this equation can be obtained using the forward and inverse Laplace transforms and 
can be represented in the form: 

a (t) = f>ie 2 ’ T erfc (— Z\ y/rj + b 2 e z * T erfc (— z 2 y/r) , (VIII-3) 
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where f>i and 62 are constants depending on initial conditions, while z 1,2 are the roots of the 
characteristic equation: 


*i,2 = -^LpV.^, ± l + W m . 


(VIII-4) 


If W' t > 0, then both roots Z\ and z<i are real and have opposite signs. Since the function 
e* 2 erfc(— z) exponentially grows for real positive z — ► 00, the steady radius a, is unstable. This 
conclusion is consistent with the result of the second-order theory (if we set e = 0 in (VIII-3)). If 
< 0 the second-order theory shows that a, is stable. In the third-order theory the stability 
criterion is W' < — \e 2 L 2 0 • This provides |argzi )2 | > \ , necessary for exponential decay of 

a (r) . In the case -\e 2 L 2 0 V?$ 20, < W, < 0, one can check using the criterion |argz li2 | > f that 

a, is stable at > 0, and a, is unstable at $20. < 0. 

The third-order approximation provides the following asymptotic expansion near the stable 

radius at large times: 


Oo(t) = a* + 


eLp V r «$2o« 
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(VTII-5) 


2. Results of Computations 

Even for a fixed substance and ambient conditions, such as for water at 1 atm ambient pressure 
and 100°C, classification of the vapor bubble behavior in acoustic fields is difficult, since three 
parameters of the pressure field u, Pa, and A p, the initial radius of the bubble introduced into the 
field, din, and the unknown 0 form a 5-dimensional parameter space. First we consider A p = 0. 
It is known from previous studies and computations [67, 38] that in this case there exist two 

equilibrium average radii of the vapor bubble, a*. 

The lower equilibrium radius is unstable and is known as the threshold radius, <Zth, since bubbles 
with initial size a{ n < ath collapse, while bubbles with Oj„ > Oth grow in the acoustic field. The 
upper equilibrium radius, a st , is stable if the stability criterion obtained in the previous section 
holds, and for bubbles with a in > a th their radii tend to a 3t at large times. This is seen from the 
phase portrait of Equation (IV-49) plotted in Figure VIII-1. The growth rate, da 0 /dt = t 2 dao/dr, 
is zero at 00 = a,. In the illustrated case the maximum growth rate is realized at small and high 
(of order 1) values of ft. 

In the case of small (3 the growth rate has a sharp peak due to the phase shift and a high 
amplitude of the resonance bubble oscillations (see Figures VII-1 and VII-2) producing non-linear 
effects. However, for (3 = 0 the rate of evaporation/condensation is zero and the growth rate is 
zero. That is why the growth rate at very small values of the accommodation coefficient, such as 
0 = 0.0001 in Figure VIII-1 is smaller than for 0 = 0.001. Such a non-monotonic dependence of 
the growth rate on the value of the accommodation coefficient was noticed earlier [2]. The case 
of high 0 demonstrates that the threshold and stable radius can be substantially smaller than the 
primary resonance radius, since they are determined by the second resonance (see Figure VII-l). 
For an acoustic frequency of 60 kHz and 0 from 0.01 to 0.1 the maximum of the growth rate is 
located between the two resonance bubble sizes. 
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It is interesting that nonequilibrium phase transitions may influence the rectified heat transfer 
at surprisingly low frequencies (kilohertz range, see Fig. VIII-2). The effect is sufficiently strong 
near the bubble resonant sizes. However, it should be studied more carefully near the second 
bubble resonance because the present weakly nonlinear theory does not cover that range and only 
the neighborhood of the primary resonance is plotted in the figure. 

Figure VIII-3 shows the relation between a* and the sound amplitude at the bubble location 
for various (3. For each [3 there exists a minimum of the plotted dependence which determines 
the amplitude threshold of vapor cavitation P For low intensity sound with Pa < 
any bubble will collapse due to the effect of surface tension. For Pa > P two equilibrium 
radii: unstable, at h (plotted by thin lines) and stable, a st (plotted by thick lines), can be found 
from Figure VIII-3. Here and everywhere below stability of the equilibrium radius was determined 
using the stability criterion in the third order approximation. However we should notice that in 
the computed cases this criterion gives practically the same results as the stability criterion in the 
second order approximation, and the transition from stable to unstable radius occurs in very close 
vicinity of the extrema of the curves. 

Figure VIII-4 shows the same dependence, but for three different frequencies. Note that 
although the present theory can be applied to high frequency ultrasonic fields it is limited by 
relatively small amplitudes (weakly non-linear approximation, e <C 1). For water at atmospheric 
pressure and A p = 0 the theory is limited to frequencies of order 100 kHz and less. For other 
substances such as cryogenic liquids the weakly non-linear approximation is valid for much higher 
frequencies. Figure VIII-5 shows that for liquid helium near the boiling point at atmospheric 
pressure, frequencies of 10 MHz and higher can be described by the present theory. Computations 
for liquid helium using the equilibrium scheme of phase transitions and the non-equilibrium scheme 
with /3 = 1 recommended at temperatures above the lambda-point [17] did not show a noticeable 
difference for frequencies up to 10 MHz. 

The curve for the equilibrium radius at 1 kHz plotted in Figure VIII-4 shows that the stable 
radius of vapor bubble for Pa ~ 0.1 atm is of order 10 cm and larger for higher Pa- This explains 
why in the numerical simulations of Hao and Prosperetti [35] performed for frequencies of about 
1 kHz, the stable radius was not achieved. The computations were carried out for millimeter 
bubbles which are much smaller than the stable size and should grow (even at very slow growth 
rates). As it is clear from Figure VIII-4 for higher frequencies stable oscillations can be reached 
for millimeter bubble sizes for 10 kHz fields and submillimeter bubbles for higher frequencies. 

Positive values of A p correspond to subcooled liquids. As seen in Figure VIII-5 for higher A p 
higher intensity acoustic fields are required for acoustic vapor cavitation. For high /? additional 
roots of the function W 2 (do) can appear, which corresponds to two unstable and two stable radii 
related to the primary and to the second resonances [38]. This is illustrated in Figure VIII-6 where 
the plotted horizontal line intersects the curve for (3=1 four times. In the illustrated case for (3 
smaller than 0.1 only one stable radius corresponding to the primary resonance is realized. 

Negative values of A p correspond to superheated liquids. In such liquids bubbles with radii 
din > d e , where a e is determined by (III-36) grow in the absence of an acoustic field. As shown in 
Figure VIII-7 acoustic fields of relatively low frequency shift the threshold of bubble growth toward 
lower sizes. In the illustrated case of 10 kHz for water at atmospheric pressure and A p = —0.01 
atm ( a e = 116 /zm) there is no equilibrium radius, and bubbles grow indefinitely if they exceed 
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the vapor cavitation threshold. It is interesting that for higher frequencies at the same conditions 
there can appear three equilibrium radii (one stable radius and two thresholds) in the range a < a e . 
In the case of 60 kHz driving frequency illustrated in Figure VIII-8, curves computed at various 
values of (3 show local minima and maxima at a < a e . 

Figure VIII-8 shows that there exists a qualitative difference between the bubble dynamics at 
low, moderate, and high (3. For example for an acoustic amplitude Pa =0.15 atm indicated by the 
dashed line at (3 = 0.04 bubbles with radii a in > 29 fim will unlimitedly grow, while for /? = 0.1 
bubbles with initial radii 27 fim < di n < 82 fim will stabilize near a & = 51 /im and only bubbles 
with a{ n > 82 fim will grow unlimitedly. The same situation takes place for (3 = 0.01, but with 
threshold radii, = 32 fim and a $ = 98 fim, and the stable radius a st = 89 fim. Such influence 
of the accommodation coefficient on the threshold and stable radii of vapor bubbles in acoustic 
fields can be exploited for measurement of /?. 

This is supported by Figures (VIII-9) - (VIII-10) which demonstrate the bubble dynamics at 
different values of the accommodation coefficient with other conditions constant. Computations of 
these curves were performed using the second order approximation for the mean bubble radius and 
the first order approximation for the amplitude of oscillation. Substantially different signatures 
of the bubble dynamics created by different amplitudes of oscillation, growth and collapse rates, 
thresholds, and stable average radii can be used for evaluation of (3. 

Figures VIII-ll-VIII-16 demonstrate the influence of initial conditions on the vapor bubble 
dynamics in acoustic fields. In all cases computations were performed using the third order ap- 
proximation for the mean radius and the first order approximation for the oscillatory part. The 
second order approximation (IV-49) provided the same results as the third order approximation 
(IV-63) where $200 — 0 was assumed. Condition $200 = 0 as well as the second order approx- 
imation corresponds to the quasi-steady initial temperature profiles in the liquid. If the bubble 
is placed into the liquid of uniform temperature several initial bubble pulsations cause increase 
or decrease of the averaged bubble wall temperature which is treated by the present theory as a 
jump at <2 = 0. Condition $200 7 ^ 0 describes this jump and can substantially influence rectified 
heat transfer. Due to the initial jump of the non-oscillatory component of temperature the mean 
bubble size starts to grow or decay proportional to \Jt. This is determined by the sign of the 
function $20 (&o) which can be negative of positive depending on Oq,P,cj, and other parameters 
(see Figure VIII-11). 

Fig. VIII-12 shows the influence of nonequilibrium phase transitions on the vapor bubble 
dynamics in a 5 kHz acoustic field. All computations were made using the numerical-analytical 
method with the same (50 fim) initial bubble radius with the only difference in (3. Computed 
average bubble radius was superposed with high-frequency oscillations obtained from the linear 
theory. It is seen that the bubble growth rates and amplitudes of oscillation strongly depend on 
(3. The maximum rates are observed near the bubble resonant radii which also are functions of /3. 

Computations presented in Figure VIIT13 and VIII- 14 show non-monotonic curves for the 
mean bubble radius for the case $200 = $20 (&tn) while for the case $200 = 0 which coincides with 
the second-order approximation computations with the same initial radius, a monotonic increase 
or decrease of the mean bubble size was obtained. It is also noteworthy that the time required 
to achieve the stable bubble size (a^ = 71 fim) in the illustrated case substantially differs for the 
two different initial temperature profiles in the liquid. However, the time required to reach the 
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state independent of initial bubble size is approximately the same for different initial temperature 
profiles (of order u}t/2n ~ 10 3 -10 4 ). We found that after reaching this state efficient computations 
for $200 7 ^ 0 can be performed using large-time asymptotic expression (VIII-5). Figures VIII-13 
and VIII-14 demonstrates that the stable bubble size can be achieved for very large times (10 7 
cycles, or about 3 minutes). However, we cannot prove this result mathematically strictly within 
the present theory (due to unknown structure of the residual terms. We show it because of a good 
agreement with straightforward computations at large times. Note that Hao & Prosperetti [35] 
reported that the bubble radius did not show stabilization at u;£/27r ~ 10 4 . 

This is an interesting effect since estimations of the characteristic time of establishment of the 
quasisteady temperature profile in the liquid are of the order t sta b ~ o,q/ki ~ 0.1 s, or 10 4 cycles. 
The explanation can be found if we repeat calculations with the initial quasisteady temperature 
profile (decaying as 1/r at infinity, which is also an assumption underlaying the second-order 
theory with respect to e). In this case the times of stabilization agree well with estimations of 
tstab, after which the average temperature exponentially reaches the steady distribution. For a 
bubble starting to oscillate in a uniform temperature field the unsteady average heat flux due to 
the initial temperature jump at the bubble surface decays as 1/2 which is determined by far 
field solution for a point heat source. Asymptotes corresponding to the principal unsteady term 
(VIII-5) are plotted in Fig. VIII-14 and they agree well with numerical results at t > t sta b- 

The temperature profiles in the liquid corresponding to the case plotted in Fig. VIII-14 (p = 
0.04) are shown in Fig. VIII-15. The dimensionless thickness of the ‘fast’ temperature sublayer, in 
the computed case was of order Lri/a 0 ~ 10 -2 . It is interesting to note that the matching boundary 
condition, $ = T 20 -b E 20 (see (IV-67)), shows a nonmonotonic behavior — first it increases and 
then decreases. This is related to a non-monotonic behavior of the mean bubble radius (see Fig. 
VIII-14) and a competition of the first (quasisteady) and the second (memory-type) terms in the 
right hand side of equation (IV-69). Fig. VIIT15 shows that even at times of order 10 6 periods 
of oscillation there exist some difference between the steady and unsteady" average temperature 
profiles. 

The possibility for the existence of steady bubble oscillations near the low-frequency resonance 
was discussed in [40 , 34] and suggested as an explanation of the observed small stable bubbles 
in liquid helium. The present theory does not predict stable bubbles at 4.2 K and atmospheric 
pressure for 52 kHz acoustic field in the range of 15 jiTri. However computations showed that 
at the parameters of the experiment the unstable threshold radius, a**, of order 15 fJ-m can be 
realized. Although this radius is unstable, the times of the instability development can reach 
several seconds, or hundreds of thousands of periods of oscillation (see Figure VIII-16; the thick 
and thin lines correspond to $200 = $20 (°in) and $ 2 oo = 0, respectively) which is of the same 
order of magnitude as the times observed in the experiments. In the experimental photographs 
some patterns of bubbles in standing waves were seen. Pattern formation in the liquid with vapoi 
bubbles is a result of strong acoustic field/bubble interaction leading to self-organization of the 
bubbles [68, 69, 70]. Such interaction as well as the bubble drift to the nodes or antinodes of the 
standing wave can also provide stabilization of bubbles of a certain size. 
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Figure VTII-1: Dependence of the growth rate of the vapor bubble on its radius for saturated 
water. The numbers near the curves show the values of the accommodation coefficient, ( 3 . The 
curve marked as “equilibrium” is computed using the quasiequilibrium scheme of phase transition. 
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Figure VIII-2: The vapor bubble mean growth rate in an acoustic field as a function of the bubble 
current mean radius at various values of the accommodation coefficient (the second order theory). 
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Figure VIII-3: Relation between the pressure amplitude and the equilibrium mean radius a, for 
saturated water. The numbers near the curves show the values of the accommodation coefficient, 
/3. The curve marked as “equilibrium” is computed using the quasi-equilibrium scheme of phase 
transition. 



Figure VIII-4: Relation between the pressure amplitude and the equilibrium mean radius a, 
for saturated water for 1 kHz, 10 kHz, and 100 kHz acoustic fields and for two values of the 
accommodation coefficient shown near the curves. 
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Figure VIII-5: Relation between the pressure amplitude and the steady radius a. for saturated 
( a 7= 0) and subcooled (Ap = 0.003 bar) liquid helium for 100 kHz, 1 MHz, and 10 MHz acousfc 

fields. 



Figure VIII-6: Relation between the pressure amplitude and the equilibrium mean radius a, for 
subcooled water. The numbers near the curves show the values of the 
0. The curve marked as “equilibrium” is computed using the 

transition. The thick lines correspond to the stable equilibrium mean radius while the thm ^ 
show the unstable (threshold) equilibrium mean radius. The dotted line [P A 0.1 ) 

each plotted curve at 2 points, except the curve for 0 = 1, which is intersected a p 
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Figure VIII-7: Relation between the pressure amplitude and the equilibrium mean radius a. for 
superheated water. The numbers near the curves show the values of the accommodation coefficient, 
(3. The curve marked as “equilibrium” is computed using the quasi-equilibrium scheme of phase 
transition. 



Figure VIII-8: Relation between the pressure amplitude and the equilibrium mean radius a, for 
superheated water. The numbers near the curves show the values of the accommodation coefficient, 
(3. The curve marked as “equilibrium” is computed using the quasi-equilibrium scheme of phase 
transition. The thick lines correspond to the stable equilibrium mean radius while the thin lines 
show the unstable (threshold) equilibrium mean radius. The dotted line [Pa = 0.15 atm) intersects 
each plotted curve at 3 points, except of the curve for /? = 0.04, which is intersected at 1 point. 
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Figure VIII-9: Dynamics of bubbles of various initial radii in superheated water at (3 = 0.01. The 
second order theory. 
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Figure VIII-10: Dynamics of bubbles of various initial radii in superheated water at ft = 1. The 
second order theory. 




Figure VIII-11: Function $ 20 (^ 0 ) describing initial conditions for the rate of the mean bubble 
radius growth in the third order approximation at various values of the accommodation coefficient 
indicated near the curves. 
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Figure VIII-12: The vapor bubble dynamics in an acoustic field at various values of the accom- 
modation coefficient (the third-order theory). The initial radius in all cases was 50 fj.m. 
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Figure VIII-13: Dynamics of the mean vapor bubble radius in water with different initial temper- 
ature profiles near the bubble, $ 200 , and different initial radii, a^, shown near the curves. The 
thick lines correspond to the temperature jump at t = 0 ($200 = $20 (&i n )) while the thin lines 
correspond to the initial quasi-steady temperature profiles near the bubble (3>200 = 0). 
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Figure VIII-14: The dynamics of the mean (period-averaged) vapor bubble radius in an acoustic 
field at various values of the accommodation coefficient and initial temperature profiles in the 
liquid. The third order theory for initial quasisteady temperature fields provide the same curves 
as the second order theory. The asymptotes show a two-term asymptotic expansion near the 
equilibrium mean radius. 
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Figure VIII- 15 : Temperature profiles in the liquid during rectified heat transfer (third-order theory, 
numerical results). The numbers near the curves correspond to the number of cycles of bubble 
oscillation. The dashed curve shows a hyperbolic profile (\I/(ao) — T^/r] for \&(ao) a t the millionth 
period of oscillation. 
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Figure VIII- 16 : Dynamics of the mean vapor bubble radius in liquid helium with different initial 
temperature profiles near the bubble, $200> and different initial radii, a in , shown near the curves. 
The thick lines correspond to the temperature jump at t = 0 ($200 = $20 i a in )) while the thin 
lines correspond to the initial quasi-steady temperature profiles near the bubble ($200 = 0 )' 
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Chapter IX. 

Analysis of Nonlinear Bubble 
Dynamics in Standing Acoustic Waves 


A two-dimensional autonomous system (V-18) is convenient for qualitative analysis. In the present 
model for fixed temperature (say water at 100°C), the degree of subcooling/superheating, A p, the 
frequency and amplitude of the acoustic field, cj and Pa, and the initial bubble radius and position, 
aoo and Xoo, can be arbitrarily varied. The value of the accommodation coefficient, /?, frequently 
is not well known for the system under consideration (since it is sensitive to trace impurities), and 
it can be considered as an unknown parameter. Finally, the value of g , can be varied over a wide 
range. These independent parameters form a 7-dimensional parametric space, which complicates 
a parametric analysis. The drag coefficient could also be varied. However, in the present analysis 
we kept it constant K p = 3, which corresponds to the drag coefficient of a spherical bubble at 
high Reynolds numbers. 


1. Zero Gravity Conditions m 

We start our analysis by consideration of the vapor bubble stability in zero-gravity conditions, 
when we can formally set g r = 0. In this case the function ^(uoj^o) has zeros at nodes and 
ahtinodes of the acoustic pressure, sin Xq cos Xq = 0. These zeros correspond to stable and unstable 
equilibrium positions of the bubble in the standing acoustic wave. The bubble is attracted to 
the pressure antinodes, xq = \ + 7 rn (n = 0,1,...), when v^{oq) > 0, and to the pressure nodes 
Xq = 7 T 7 i, when ^ 22 (^ 0 ) < 0. We will call the bubble ‘A-stable’ if it is attracted to the pressure 
antinode and ‘N-stable* if it is attracted to the pressure node. We also define the bubble as V- 
stable’ if its average volume (or a$) is stable in an acoustic field. It is clear that for overall bubble 
stability both positional and volume stability should be achieved. Note that the N-stable bubble 
will grow or collapse since the local acoustic pressure is zero at the pressure node and the bubble 
is V-unstable in the absence of the acoustic field. This is consistent with Eq.(V-19) that gives 
W 2 (ao,7rn) = ^ 20 (^ 0 ), where w 20 M is a monotonic function. Therefore, only A-stable bubbles 
corresponding to ^ + 7rn) = ^ 20 (^ 0 ) d" ^22 ( a o) can he V-stable and generally stable. 

The function ^ 22 (^ 0 ) does not depend on A p and P^. Thus, at fixed ambient conditions, 
and known drag coefficient, K AZ , the zeros of this function, a po3) depend only on the acoustic 
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frequency and the accommodation coefficient, Op 03 = ap OS (uj,P). The function W^Oo, \ -fi 7 m) 
depends on the same parameters as well as on A p and Pa . For the zeros of this function we have 
a vo i = a vo i (u>, (3, A p, Pa) ■ Consider the following three cases showing qualitatively different results 
about the roots a vol . 


a. Saturated Liquids 

In this case the ambient pressure corresponds to the saturation line, and A p = 0. Two roots 
of equation ^(ooi f + 7Tn ) ~ 0 appear only if the amplitude of the acoustic field exceeds the 
threshold of vapor cavitation Pa > Pa* ( w , /?) (Figure IX-1). The value of Pa* (w, (3) can be found 
by excluding oq from the simultaneous conditions: 


W 2 {ao, | + 7m; P A , u, fi) = 0, 


dP A (ao,w,0) _ 0 

da 0 


(IX-1) 


Curves Oo = a vo i (u, (3, 0, Pa* (tu, /?)) determine a region where stable bubble volume oscillation 
can exist. Projection of this regions onto the (a, u) plane for water for the value (3 — 0.04 is shown 
in Figure IX-2a. It is interesting to note that the regions of A-stability and V-stability in the 
computed cases do not intersect. In other words, for A p = 0, water vapor bubbles at 100° C are 
unstable in any acoustic field. Unfortunately we cannot prove this observation for other substances 
or for water in other conditions because of the high dimension of the parametric space. 

Figure IX-3 demonstrates phase trajectories of the dynamical system (V-18) computed at Pa > 
Pa* (oj,/3) . The phase portrait is periodic with the period equal to half the acoustic wavelength. 
The figure plots only one period. All trajectories end at Oo = 0, despite the existence of regions 
of bubble growth. However, all growing bubbles finally move toward the pressure node, where the 
acoustic field cannot support the bubble growth and the bubble collapses. Note that only bubbles 
exceeding a critical size can grow near the pressure antinode. 

There exists a possibility for bubble average radius oscillations (spiral trajectories). The mech- 
anism of such oscillations is the following. A bubble growing near the pressure antinode is forced 
to move towards the pressure node. There, the local amplitude of acoustic pressure is lower and 
it collapses. At smaller sizes the direction of the acoustic radiation force changes and the bubble 
moves back toward the pressure antinode where it can grow again and repeat the cycle. However, 
our analysis shows that these oscillations are unstable and the bubble will always collapse. Note 
that such oscillations can occur in both superheated and subcooled liquids. An example of the 
oscillating average bubble size behavior is shown in Figure IX-6. 


b. Superheated Liquids 

This case corresponds to A p < 0. Any bubble exceeding the critical size a e , Eq.(III-36), will grow 
even in the absence of the acoustic field. In contrast to the previous case, now three roots of 
equation Wj(oo, f + 7 m) = 0 can exist (see Figure IX-lb). This determines two branches of the 
curve a 0 = a vol (u>, /?, A p, Pa* (w, {3, Ap)) bounding the region of V-stability. Threshold values for 
Pa* A p) are determined by Eq.(IX-l). Projection of this region onto (a, u) plane for water 

at A p = —0.01 atm and (3 = 0.04 is shown in Figure IX-2b. As in the previous case the regions 
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of A-stability and V-stability for water vapor bubbles at 100°C do not intersect and, therefore, 
bubbles are unstable in superheated liquids in any acoustic field. 

Figure IX-4 displays the phase trajectories of the dynamical system (V-18). All trajectories 
either end at do = 0, or grow unlimitedly corresponding to bubble growth due to the liquid 
superheat. Figure IX-7 illustrates the dynamics of two bubbles with the same initial position in 
the wave, but slightly different initial sizes near the critical size a e . Both of them move relatively 
quickly towards the pressure antinode (the x-coordinates of the centers in these cases practically 
coincide). At the antinode the smaller bubble collapses, while the larger bubble grows. This is 
determined by the bubble size with respect to a e at the moment of reaching the pressure antinode 
of the standing wave. 

c. Subcooled Liquids 

In this case A p > 0. The intriguing feature of this case is the possibility of existence of four 
roots of the equation W 2 {a 0 ,f + tto) = 0 (Figure IX-la). This determines three brandies of 
curve oo = (w, P, A p, P A . (v, P, A p)) bounding two regions of V-stability. Projections of these 

regions onto the {a,u) plane for water at A p = 0.02 atm and P = 0.04 and P = 1 are shown in 
Figure IX-2c. It is seen that an intersection of the regions of stability exists, which provides both 
A- and V-stability. Therefore, in subcooled liquids vapor bubbles can be stabilized by an acoustic 
field. However, we should notice that the amplitude of this field is constrained by two similar 
values of P At corresponding to the maximum and minimum of the function P A (a 0 ). For example, 
in the case shown in Figure IX-la at A p = 0.02 atm, the bubble can be stabilized in a narrow 
range of amplitudes 0.0885 atm < P A < 0.0892 atm. The range can be wider at different ambient 
conditions or for substances other than water. 

The region of stability occurs at relatively low frequencies, when the primary resonance radius 
and the condensation-evaporation (or ‘second’; Marston, 1979; Hao and Prosperetti, 1999) reso- 
nance radius are substantially separated. The stabilization occurs near the second resonance. Note 
that Marston and Greene (1978) observed stable bubble oscillations in subcooled liquid helium 
and explained them as an effect of the second resonance. ^ - 

Figure IX-5 shows phase trajectories of the dynamical system (V-18) in the neighborhood of 
the stable state. All trajectories from the attraction region end at the point corresponding to the 
stable state. Note that this region is limited and is relatively wide for bubble size and relatively 
narrow for bubble position. Therefore, to obtain stable oscillations the bubble should be initially 
situated close to the pressure antinode of the acoustic wave. The process of stabilization of bubble 
size and position is shown in Figure IX-8 for two bubbles with initial sizes smaller and larger than 
the stable radius and the same initial position. An interesting feature is that the bubble of the 
larger initial size contracts below the stable radius and then expands. This is consistent with t e 
phase trajectories of Figure IX-5. 


2. Effects of Gravity 

According to Eq.(V-18), gravity influences the bubble growth rate and the average translational 
velocity. However, the effect of gravity on bubble growth rate is rather small (if it is taken into 
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Radius (pm) Radius (nvn) 

Figure IX-1: The amplitudes of 1 kHz and 100 kHz isotropic acoustic fields at which a bubble can 
experience steady (stable or unstable) oscillations for water at /? = 0.04 and T = 373 K. 

account just in the last term of Eq.(III-4)). The most important effect is destabilization of the 
bubble position. 

As an example of destabilization, Figures IX-9-IX-10 show phase trajectories of system (V-18) 
at various levels of gravity. In zero gravity conditions the system has a stable state (Figure IX-5). 
If some level of gravity is imposed the equilibrium position shifts from the pressure antinode. The 
larger gravity is, the larger is this shift. This takes place until some threshold value of gravity 
corresponding to reaching the boundary of the attracting region. For gravity higher than this 
threshold, the steady state solution longer exist. Figure IX-11 demonstrates the bubble dynamics 
in subcooled liquids. 

It is interesting that gravity can turn unstable oscillations of average bubble size and position 
into stable oscillations. In zero gravity the system has unstable equilibrium points (see Figure IX- 
2a). For increasing magnitude of gravity (Figures IX-12-IX-14), the two illustrated equilibrium 
points are not symmetric, and a limit cycle appears near the left equilibrium point. However, 
further increase of gravity up to one g destabilizes this limit cycle. Figure IX-15 demonstrates the 
bubble dynamics in saturated liquids. The stable oscillations of the average bubble radius and 
coordinate of the center can also be observed in subcooled and superheated liquids. 

Figures IX-16-IX-17 present computaions of the mean bubble radius and position for different 
values of the accommodation coefficient and different levels of gravity for all other parameters 
fixed. It is seen that the measurements of the accommodation coefficient can be performed at 
zero gravity, reduced gravity, and normal gravity conditions. However, for each level of gravity 
there exist specific patterns or regimes of bubble dynamics which should be taken into account in 
selection of optimal conditions for measurement of the accommodation coefficient. 




Frequency. kKr Frequency. kHr 

Figure IX-2: The regions of A- and V-stability of bubbles at zero gravity for saturated, super- 
heated, and subcooled water at /? = 0.04 and T = 373 K. 



Figure IX-3: The phase trajectories for saturated water at Og, /? = 0.04 and T — 373 K. 10 kHz 
acoustic field, € = 0.099. 



Figure IX-4: The phase trajectories for superheated water at Og, Ap = -0.01 atm, (3 — 0.04 and 
T = 373 K. 60 kHz acoustic field, e = 0.123. 
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Figure IX-5: The phase trajectories for subcooled water at Og, Ap = 0.02 atm, /? = 0.04 and T 
373 K. 1 kHz acoustic field, e = 0.088. 
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Figure IX-6: The bubble dynamics corresponding to a phase trajectory from Figure 3. 
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Figure IX-7: The bubble dynamics corresponding to phase trajectories from Figure 4. 
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Figure IX-8: The bubble dynamics corresponding to phase trajectories from Figure 5. 



Figure IX-9: The phase trajectories for subcooled water at 0.05g, Ap — 0.02 atm, (3 — 0.04 and 
T = 373 K. 1 kHz acoustic field, € = 0.088. 



Figure IX-10: The phase trajectories for subcooled water at Ig, Ap = 0.02 atm, /? — 0.04 and T 
= 373 K. 1 kHz acoustic field, e = 0.088. 


97 



Figure IX-11: The bubble dynamics in subcooled water at various levels of gravity Ap = 0.02 
atm, /? = 0.04 and T = 373 K. 1 kHz acoustic field , e = 0.088. 



Figure IX-12: The phase trajectories for saturated water at 0.05g, (3 = 0.04 and T = 373 K. 10 
kHz acoustic field, e = 0.099. 



Figure IX-13: The phase trajectories for saturated water at O.lg, (3 = 0.04 and T = 373 K. 10 
kHz acoustic field, e = 0.099. 
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Figure IX-14: The phase trajectories for saturated water at lg, P = 0.04 and T — 373 K. 10 kHz 
acoustic field, e = 0.099. 



Time, s 


Figure IX-15: The bubble dynamics in saturated water at various levels of gravity, P = 0.04 and 
T = 373 K. 1 kHz acoustic field, e = 0.088. 
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Figure IX-16: Dynamics of a vapor bubble in a standing acoustic wave at different values of the 
accommodation coefficient and normal gravity conditions. 



Figure IX-17: Dynamics of a vapor bubble in a standing acoustic wave at different values of the 
accommodation coefficient and reduced gravity conditions. 
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Figure IX- 18 : Dynamics of a vapor bubble in a standing acoustic wave at different values of the 
accommodation coefficient and zero gravity conditions. 
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Chapter X. 
Conclusions 


Our study conducted using analytical and numerical methods shows that the kinetics of phase 
transitions can affect the dynamics of vapor and vapor-gas bubbles over a broad range of frequen- 
cies and bubble sizes. For water at atmospheric pressure, the frequency range 1-100 kHz and the 
bubble size range 10 lira - 10 mm were investigated. During the study several important physical 
effects were uncovered, including 

• Strong influence of the kinetics . of phase transitions on the vapor temperature, 

• Influence of the kinetics of phase transitions on the dynamics of vapor-gas bubbles with 
substantial inert gas content, 

• Influence of the Soret-Dufour effect on bubble resonance oscillations, 

• Memory-type thermal effects in the liquid during rectified heat transfer, 

• Effects of slow bubble mean radius/position oscillation and stabilization in standing wave. 

These effects can be exploited for proper design of experiments. Such experiments would 
be performed prior to finalizing the design of an instrument to investigate practical issues. We 
found that there exist optimal superheats, concentrations, frequencies, and mean bubble sizes 
which provide the highest sensitivity of the bubble dynamics to the value of the accommodation 
coefficient. This sensitivity can be substantially reduced in conditions that are far from optimal. 
The experimental setup should allow for variation and control of the liquid temperature and the 
pressure in the cell. It is also important to select appropriate methods of bubble generation and 
measurement of bubble size and position. Depending on the methods employed the optimal ranges 
of bubble sizes and acoustic frequencies can be determined for measurement of the accommodation 
coefficient. The characteristic acoustic frequency will determine the size of the setup (Helmholtz 
resonator for standing waves), which is also should influenced by the methods of bubble size 
and position measurement and bubble generation. The experimental design also account for 
variation in the optimal sensitivity of the bubble dynamics with the order of magnitude of the 
accommodation coefficient. If is to be measured over a broad range of expected values (say 
0.001-1 for water), then this range can be subdivided into several bands (say 0.001-0.01, 0.01-0.1, 
and 0.1-1). For each band, an optimal set of conditions can be found. 
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